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Abstract 



Present work is a theoretical study on the stability of the thermotropic biaxial nematic liq- 
uid crystal phase in model systems. Its main aim is to present the phase diagrams of spatially 
uniform liquid mesophases and to identify the molecular parameters that influence the stability 
of the biaxial nematic. The diagrams are obtained by means of the Local Density Functional 
Theory in low density approximation, and the relation between the molecular parameters of 
the models and macroscopic properties of the system close to the transition point are obtained 
by means of bifurcation analysis. We consider three model systems; the so-called L = 2 
model (the lowest coupling model of the orientational part of pair potential), the biaxial Gay- 
Berne interaction, and the bent-core system. For the second one, we also briefly investigate 
the elastic constants and comment on the smectic phases. In every case we take into account 
rigid molecules. 

Firstly, we study a generalized version of Straley model, which should be considered as 
the simplest one giving rise to the biaxial nematic phase; it is to the biaxial nematic what 
Maier-Saupe model is to the uniaxial nematic. Its analysis allows to construct a generic, ex- 
act in mean field bifurcation diagram as function of potential parameters. By considering the 
symmetries of the L = 2 expansion of pair interaction, valid also for pair direct correlation 
function, we find the duality transformation which connects the states at different temper- 
atures. The so-called self-dual points, i.e., the points in space of potential parameters left 
invariant under the action of the duality transformation, coincide with Landau points. The 
following Landau region is found as well as the tricritical points. Despite its simplicity, the 
model exhibits a rich behaviour. Many systems can be accurately approximated using this 
approach, and located on the presented diagram. 

The second model studied is a Gay-Berne potential generalized to the soft ellipsoids of 
three different axes. It showed the biaxial nematic in Monte Carlo study, and therefore it gives 
us the chance to confront our results with simulations and study the interaction parameter 
space in more detail. We also compare the bifurcation diagram following from Local Density 



Functional Theory with transition points acquired by the minimisation of the Helmholtz free 
energy. In comparison with the simulations, the approach used overestimates the transition 
temperatures. We study the interaction further and locate the Landau points (also called self- 
dual) induced by the increase of both the shape and energy biaxiality. In the former case, 
our results indicate that the introduction of attractive forces slightly shifts the position of the 
self-dual point from its location for hard biaxial ellipsoids. We find that the direct isotropic - 
biaxial nematic transition occurs at Landau points, and show how the concurring biaxialities 
influence their position. The results indicate that for this model self-dual region predicted for 
hard, biaxial ellipsoids by the so-called square root rule retains its significance, and provides 
qualitatively correct estimations of Landau points positions. 

Analysis of the bent-core system is aimed at studying the shape induced effects and the role 
of the dipole-dipole interaction. The model molecules are built using two and three Gay-Berne 
interacting, prolate, soft ellipsoids of uniaxial and biaxial symmetry. We study the influence 
of dipole-dipole interaction, in case of two and three uniaxial arms, by introducing a trans- 
verse dipole moment lying along the C 2 symmetry axis, in the plane of an angle between the 
arms (opening angle). For non-polar case of two uniaxial arms, we find that by changing the 
opening angle we can obtain a Landau point at 107°, in agreement with previous studies for 
hard bent-core molecules. Surprisingly, for this model the inclusion of attractive forces does 
not influence its position. Once the arms deviate from the uniaxial symmetry and become 
biaxial ellipsoids interacting via biaxial Gay-Berne potential, the self-dual point evolves to the 
line of Landau points ranging in the opening angle between 121° and 128°. In both cases, the 
opening angles at which self-dual points occur, differ from the experimental estimates of 140°. 
However, a study on system of bent-core molecules with opening angle of 90° was published 
recently, and those results remain in qualitative agreement with our analysis of the model of 
three uniaxial Gay-Berne parts where we find the bifurcation diagram with the isolated Landau 
point at the angle between the arms equal to 89°. The introduction of weak dipoles (relative to 
the strength of the Gay-Berne energy) leads to the shifting of the Landau point towards lower 
angles. For three uniaxial parts model, moderate dipole magnitudes provide the diagram with 
a line of self-dual points, range of which in opening angle increases with dipole strength. The 
maximum length of the interval of the opening angle for which a direct isotropic - biaxial 
nematic transition occurs is 23° between 63° and 86°. For the strongest dipoles studied, that 
line begins to shrink and shift towards higher angles. In bent-core models, again, the Landau 
points are the only places where isotropic - biaxial nematic transition takes place, the novel 
feature is the appearance of the self-dual line. 

Since the elastic constants are important quantities in the theory and applications of 
liquid mesophases, in final chapter we present the preliminary studies on temperature depen- 



dence of a set of bulk biaxial elastic constants for the biaxial Gay-Berne interaction in L = 2 
model. Our results are in agreement with previous findings for the lattice biaxial model; we 
recover the splay-bend degeneracy, and find that the constants associated with one of three di- 
rectors in the biaxial nematic are always negative. Also, in the rod-like molecular regime the 
relative values of the uniaxial elastic constants are higher than biaxial ones. This behaviour 
changes in vicinity of the Landau point, where the constants associated with one of the biaxial 
directors become dominant. 

Finally, we make a comment on the smectic-A phases. Although the thesis is devoted to 
the spatially uniform liquids, the competition between the smectic order and biaxial nematic 
should be considered since the former may gain stability earlier and prevent the formation of 
the spatially uniform biaxial ordering. In order to address this issue, and as an example, we 
limit ourselves to the orthogonal smectic-A phases. We present the behaviour of a complete set 
of leading order parameters for the biaxial Gay-Berne potential, based on the self-consistent 
equation for equilibrium one-particle distribution function in L = 2 model, including uniaxial 
and biaxial smectic-A states. The temperature dependence of order parameters is presented 
for two sets of interaction parameters. One of them, namely the one where molecules are 
strongest attracted to their sides (the lateral forces are strongest), was the only case for which 
Monte Carlo simulations discovered a stable biaxial nematic. We find that in this case the 
temperature range of the lower symmetric nematic state is significantly wider than for the 
model of strong face-to-face attractions. For fixed density and with increasing temperature, 
biaxial smectic is found to melt to biaxial nematic which in turn transforms to uniaxial ne- 
matic and for higher temperatures the system undergoes a phase transition to isotropic liquid. 
This phase sequence is in agreement with simulations results. In final section, we extend the 
approach used for spatially uniform states to include the orthogonal smectic phases of uniaxial 
and biaxial symmetry, and we present the additional bifurcation equations. 
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Chapter 1 
Introduction 



A hundred and twenty years ago it became clear that between the usual fluid and solid crystal 
there exists an intermediate phase of matter. Substance in the new phase remained liquid, 
however, system exhibited extraordinary anisotropic optical properties, which could be ex- 
plained by microscopic studies. It was found that in the new states a degree of order was 
formed on the molecular level. Apparently the molecules no longer possessed random ori- 
entations, as in the normal liquids, but were on average oriented along one, system-wide 
direction. A certain subclass of those intermediate phases is a subject of the present work 
for which this chapter serves as an introduction. We will begin with a brief description of the 
historical background and review of the research conduced so far, then continue to introduce 
the matter of the thesis, its aim, and the employed methods. 



1.1 Discovery of ordered liquids 



The first indications of the new phase of matter are due to W. Heintz who in 1850 studied 
stearin. He found that the substance apparently possessed two melting points, and that was a 
new phenomenon altogether. With the increase of temperature, it firstly became cloudy, then 
opaque, and eventually turned into a clear liquid. Later, in 1854, Rudolf Virchow described a 
soft, floating substance present in the human body, protecting the nerve fibres, which he called 
myelin [1] (see Fig. 1.1). By the work of Carl Mettenheimer three years later this substance 
was found to show birefringence. Those were the first indications of the new phenomena. 
Today the systems studied by Virchow and Mettenheimer would be called lyotropic, since 
concentration of water is a driving force for the phase transitions. In contrast, the materials for 
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which this role is played by temperature would be called thermotropic 1 . Surprisingly, it took 
over thirty years to realize the meaning of those findings. It is the year 1888 that is recognised 
as the year of the discovery. At that time Friedrich Reinitzer synthesized cholesteryl benzoate 
(C 3i H 50 O 2 , see Fig. 1.2) and observed with his microscope iridescent colours as well as simi- 
lar temperature driven behaviour as Heintz did 38 years earlier. Intrigued by the phenomenon 
he consulted Hofrath von Zepharovich, and by his advice wrote a letter to Otto Lehmann on 
14th of March, 1888 [2]. Otto Lehmann, born a year before Virchow's discoveries, at the 
time of receiving the letter was a young associate professor, although already a well known 
figure in the field of phase transitions phenomena. He just finished his "Molekularphysik", an 
impressive work of two volumes enclosing more than 1500 pages with 620 figures and over 
2000 citations. Friedrich Reinitzer, born 1857, at the time was a botanist. His conscientious 
observations [3] and correct conclusions were one of the main factors that helped to realize 
the meaning of the discovery. He was the right man in the right time and place, and a lot of 
credit for the early discoveries is rightfully attributed to him. On March 14th 1888 he wrote: 

The substance has two melting points, if it can be expressed in such a manner. 
At 145. 5°C it melts to a cloudy, but fully liquid melt which at 178.5°C suddenly 
becomes completely clear. On cooling a violet and blue colour phenomenon ap- 
pears, which then quickly disappears leaving the substance cloudy but still liquid. 
On further cooling the violet and blue colouration appears again and immediately 
afterwards the substance solidifies to a white, crystalline mass. 




Figure 1.1: Drawing of myelin as shown by O. Lehmann, after [4]. 

Those intermediate states, ordered and yet liquid, at first were called by Lehmann "Fliessende 
Kristalle" that is "flowing crystals", which later evolved to "liquid crystals" or mesoph&ses 

'The distinction is a bit more subtle, since in lyotropic materials the temperature is also important, but defi- 
nitely for thermotropic phases the concentration of water is not the issue. 
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from a Greek word "meso" meaning "middle" or "intermediate". It was the cooperation of 
Reinitzer and Lehmann that marked the first chapter of the liquid crystals history, which is 
recognized to began on 14th March 1888, the day the letter was sent which opened a new field 
of physics. 

Today there is a great amount of known liquid crystalline structures, only a brief descrip- 
tion of those would consume many pages, if not a bundle of volumes. The mesogenic sub- 
stances play an important role in our world. They are not only essential for a widely available 
and most useful LC displays, but also, as was mentioned, were found early on in living organ- 
isms; besides myelin, the living cell membranes manifest the behaviour of the lyotropic liq- 
uid crystals. The overall experimental issues, technological applications and interesting chal- 
lenges in theoretical description make those materials an important part of modern physics. 

The reasons for the emerging of liquid crystalline ordering are one of the fundamental 
problems that have been addressed many times. There exists a magnificent variety of sub- 
stances that manifest liquid crystalline behaviour, and yet many more can be synthesized, and 
there are numerous factors that can be responsible for mesophase formation. At least in most 
cases some anisotropy of shape of molecules is a required condition. Another is the anisotropy 
of the intermolecular forces. Those two basic issues can be set in the centre of the problem 
(together with others of which we say nothing). 

One of the interesting problems is the question of the conditions which are required for 
the stabilization of a liquid crystal structure, like temperature, density, pressure, or external 
fields. We can assume that we have at our disposal the microscopic quantities like molecular 
shape and characterisations of forces, and we want to find the macroscopic, equilibrium prop- 
erties of the system. The perfect way to establish the relation between those is the statistical 
description. This work is an example of such approach. 




Figure 1.2: Cholesteryl benzoate. The first synthesized liquid crystalline substance. 
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1.2 Nematic order 



Most of the mesogenic substances consist of molecules that do not possess any definite 
shape. However, in a given liquid crystal state the molecules are rotating much faster then they 
are moving through the sample, i.e., the time scale of the rotation is few orders of magnitude 
different from the time scale in which a measurement takes place. In effect the molecule acts 
as an object of certain averaged shape. Therefore we can approximate the complicated struc- 
ture of the compound by a simpler object of a given symmetry and think of the molecules as if 
they possessed a defined shape. In particular, many liquid crystals are built from molecules of 
average shape resembling a rod, which can be approximated by a prolate ellipsoid of revolu- 
tion. In the system of elongated, rod-like molecules in high temperature or low density regime 
we will observe the usual, disordered fluid which is traditionally called an isotropic phase 




Figure 1.4: Snapshot of molecular structure of disc-like (on the left) and rod-like 
nematic (on the right). Director n is marked. 
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(I so). By lowering the temperature, or increasing the density, it can evolve to a more ordered 
structure, possibly a liquid crystalline phase. The simplest of those is called a nematic phase, 
from a Greek word "nemato" meaning "thread". The complete disorder of isotropic liquid in 
the nematic state is replaced by a tendency of long molecular axes to on average align along 
given direction, called a director. The long range correlation of orientational degrees of free- 
dom appears while the full translational symmetry is upheld. Because there is only one system 
wide, distinguished direction, the state is referred to as uniaxial nematic (N^). Snapshots of 
the system exhibiting the isotropic and nematic states are depicted in Fig. 1.3. This state is 
additionally characterized by the invariance with respect to the rotation around the director, 
and reflection in the perpendicular plane, and therefore it is said to be of AxA symmetry, since 
this group contains those transformations. In present work a phase or molecule is considered 
to be uniaxial if it is left invariant under the operation of D^h symmetry group. 

The model molecular uniaxial symmetry can be realised in another way. Some mesogens 
consist of molecules that can be well approximated by oblate ellipsoids of revolution, resem- 
bling a disc. The uniaxial nematic phase in a system of those objects is characterized by order 
of shorter axes. It is called a disc-like or oblate nematic (N[/_), in opposition to formerly 
described rod-like or prolate nematic (N{/ + ). Both phases are shown in Fig. 1.4. Naturally, 
these states may be realized in a systems of molecules with more complex shapes, inclusion 
of which may bring a possible new behaviour of the system. 

We can consider a deviation from the D^h symmetry and take into account the model 
molecules possessing the shape of ellipsoid with three different axes. If we choose one of 
them to be significantly longer, relative to the remaining two, we can obtain the N u+ resulting 
from the ordering of the longer axes. While the system is in uniaxial nematic phase, we can 
wonder if it is possible to introduce the correlations of shorter axes, while maintaining a lack 
of long-range order of molecular centres of masses. In that way we would acquire a second 
director in the plane perpendicular to the uniaxial one, and third perpendicular to those two. 
The symmetry of this state consists of reflections in three perpendicular planes and is denoted 
by D 2 h ■ This phase is called a biaxial nematic (N B ) and is represented in Fig. 1.5. As of 
now, it is a "hot subject" in the soft matter field from both theoretical and experimental point 
of view. On the experimental side, it is still unclear why some materials, especially bent-core 
and tetrapode systems, exhibit spatially uniform biaxiality while the others do not, or what 
molecular parameters are required to stabilize the biaxial nematic. In the light of recent dis- 
coveries, it seems that we are getting closer to obtaining the answer to these issues, and once 
they are resolved, the technological applications will be most promising. In this respect the 
theory can be closely related to experiment, since it can point out the way which should be 
undertaken. is one of the two simplest mesophases, yet the introduction of two directors 
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in addition to the uniaxial one brings a rich variety of new behaviour. Present work is devoted 
to that state. In the following sections we present a description of the research conduced so far 
in context of the biaxial nematic phase. 



1.3 Biaxial nematics 



1.3.1 First indications of stable biaxial nematic phase 

The lower symmetric nematic phase was predicted for the first time by Freiser [5] in 1970. 
Essentially he generalized the Maier-Saupe model [6] to non-cylindrical molecules. In the 
mean field approach, he found the first order transition from isotropic liquid to the uniaxial 
nematic phase to be followed by a second order transition to the biaxial nematic. 

One year later another approach to the biaxial phase was presented. The Landau [7] de- 
scription of phase transitions was extended by de Gennes and shown [8] to produce a stable 
biaxial nematic. 








1 





Figure 1.5: Snapshot of molecular structure of biaxial nematic (on the right) with two 
directors marked as arrows, the third one is perpendicular to the surface of the pic- 
ture. On the left an example of the biaxial molecule is shown and three perpendicular 
planes; reflections in those constitute the D 2 h symmetry. 



In 1973, Alben proposed a different method of stabilisation of the biaxial ordering [9]. It 
was based on an assumption that in a system of molecules of oblate and prolate shape the two 
types of nematic ordering and N[/ + can form a mixture, i.e., the system would not exhibit 
any areas rich in rods and with low concentration of disks or vice-versa. Furthermore, the two 
directors should then be perpendicular and in effect biaxiality could emerge (see Fig. 1.6(a)). 
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Alben studied the lattice model of hard molecules, where only steric interactions are taken 
into account. The biaxial phase was found to be stable, and at a certain concentration of discs 
even a direct I so - Ng transition took place. An example of the diagram from that study is 
presented in Fig. 1.6(b). 

A year after Alben's work, another study was presented; following from the success of 
Maier-Saupe model for isotropic - uniaxial nematic transition Straley proposed a generaliza- 
tion of the orientational part of their potential to the case of objects of D 2 h symmetry [10]. 
The new two-point interaction had four parameters, which were related to the dimensions of 
rectangular blocks by requiring the consistency with excluded volume present in Onsager's 
theory. The potential parametrized in this way was then used to construct a mean field the- 
ory. The resulting phase diagram revealed a first order Iso - Nu transition which in lower 
temperatures was followed by a second order Nu - N# transition. Also it exhibited a direct 
transition between isotropic liquid and biaxial nematic at self-dual point, as can bee seen in 
Fig. 1.7. In this work Straley also for the first time presented a complete set of four molecular 
order parameters needed to describe uniaxial and biaxial nematics. 

Most of the above theoretical studies in comparison to real systems can be considered 
huge simplifications. Not only many microscopic aspects of liquid crystal phenomena, like 
intermolecular forces, were neglected or crudely approximated, but also some structures, like 





100 



(a) Picture of mixture of rods and discs. If 
the two uniaxial phases Ny + and Nj/_ mix 
two directors fii and can be perpendic- 
ular. 



(b) One of the phase diagrams presented in 
[9]. Effective temperature T* as in [9] ver- 
sus discs concentration for elongation of pro- 
late molecules equal to 5 and width of oblate 
ones 3. 



Figure 1.6: Schematic picture of rods and discs giving rise to biaxial nematic phase, 
and example of phase diagram for model mixture of rods and discs (on the right). 
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Figure 1.7: Phase diagram adapted from [10], breadth B as function of dimensionless 
temperature t. Length and width of rectangular blocks are set to 10 and 1 respectively. 
I so - N B transition occurs at self-dual point B = y/W. 

smectics, were not included. Nevertheless, those models provided a stable biaxial nematic 
and predicted no particular difficulties in the stabilisation. Since the fundamental problems 
in the formation of biaxial ordering, had they existed, should have emerged in those studies, 
it came as a surprise that the biaxial nematic was much more difficult to observe than other 
mesophases. It took more than ten years from the first prediction to the first experimental 
realisation of Ng. In 1980, Yu and Saupe published a report on the spatially uniform biaxial 
ordering in lyotropic ternary system of potassium laurate-l-decanol-D 2 [11]. By means 
of microscopic studies and deuteron resonance (NMR), two uniaxial nematic states denoted 
NlO^u-) an d A r c(N;7 + ), and the biaxial phase were found. Phase diagram from that study is 
presented in Fig. 1.8. For D 2 concentration above 68 wt. % the system exhibits only the N L 
phase, characterized by micelles of disc-like shape. When D 2 concentration is lowered to 
67.8 wt. % the phase Nq forms, with rod-like micelles. Between those two states, the biaxial 
nematic phase occurs which on heating and on cooling is transformed into the N L state. 

1.3.2 Further theoretical studies 

Mixtures 

Study by Alben work rendered the binary mixtures of rod-like and disc-like molecules as 
a possible candidate for the realization of a stable biaxial nematic phase in thermotropic sys- 
tems [9]. In the early 80s, this matter was re-examined. Models of mixtures of hard molecules 
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Figure 1.8: Phase diagram of potassium laurate-l-decanol-D 2 as presented in [11]. 
Uniaxial nematic N L at a certain concentration of D 2 transforms to the biaxial ne- 
matic phase N B . 



of equal volume were shown to produce a stable biaxial phase [12, 13]. The effects of ex- 
cluded volume which were studied in Onsager's approximation lead to a first order transition 
from isotropic to uniaxial nematic phase, and to a second order Njy - N B transition [12]. This 
work also confirmed that point of a direct I so - Nb transition can be found. A similar study 
on the extension of Onsager theory showed as well a stable biaxial nematic phase at a certain 
concentration of molecules and fraction of discs and rods [13]. The inclusion of both long 
and short range interactions in a mean field approach using a type of Van der Waals theory 
also gave Ng[14]. In this study the system was divided into cells of a given volume; within a 
cell molecules adopted three discrete orientations, and interacted via a short-range repulsive 
potential, while the long range attractive forces acted between the cells. 

The predictions of the above studies [12, 13, 14], as well as of Alben [9], are correct, 
provided the two uniaxial nematic phases of discs and rods mix. It is not obvious however, if 
the system does not undergo a demixing transition before the formulation of the biaxial ne- 
matic phase. The configuration where states N[/_ and N u+ axe separated, i.e., in the system 
appear regions rich in molecules of one type, may gain stability instead of N B . This issue was 
addressed by Palffy-Muhoray and de Bruyn [15]. They applied previously developed [16], 
generalized Maier - Saupe model in mean field approximation to binary mixture of rods and 
discs of equal volume and found it phase separates when the arithmetic mean rule was used 
for interactions between different types of molecules [15]. Similar model was developed to 
study the effects of external fields and induced biaxiality [17]. Again, the system exhibited 
demixing and the biaxial nematic phase lost stability at the expense of two separated uniaxial 
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nematics. 

Studies mentioned above turned the attention to the question of the required conditions 
for biaxial phase to become stable against mixture separation. For the Palffy-Muhoray and 
de Bruyn model [15] it turned out that in order to stabilise the biaxial nematic phase against 
demixing it was enough to remove the constraint of the arithmetic mean for the inter-type 
potential strengths. It was also noted that the demixing phenomena was related to the interac- 
tion strengths between unlike molecules [18]. This issue was partly addressed in [19] where 
Onsager-like theory (variational cluster approximation) was used to identify the molecular 
size, elongation, and interaction strength ratio between different kinds of molecules that make 
biaxial nematic phase possible. There it was suggested that the increase in rod-disc interac- 
tions can be of importance and should be taken into account. Later, on grounds of Onsager 
theory, which was generalized in order to include selective attractions between rods and discs, 
it was demonstrated that these interactions could stabilize the Ng against demixing [20]. Sim- 
ilar results were presented for a mixture of biaxial, D 2 h - symmetric molecules in [21], where 
within mean field approximation it was proven that with high enough biaxiality and inter- 
action strength Ng was stable. Changing the relative potential strengths between the unlike 
molecules brings a degree of asymmetry to the mixture. In terms of Onsager approximation 
the asymmetric systems of oblate and prolate molecules can be studied by changing the ratio 
between excluded volumes of rods and discs. Such an attempt was made and several demixing 
scenarios were demonstrated, and at a certain excluded volume ratio the biaxial nematic was 
proven to be stable [22]. 

In Onsager model for I so - transition the orientational entropy, which is maximal in 
the isotropic phase, competes with the entropy of packing. At high enough density the latter 
can "win" and the global minimum of the free energy can be associated with the orientation- 
ally ordered state. In the case of binary mixtures another term is added, which corresponds 
to the entropy of mixing, which is maximal when the two phases are mixed. The competition 
of these three terms may lead to stabilisation of the biaxial nematic phase. However, when 
the two components do not mix it is possible for the entropy of mixing to "loose" and for the 
system to separate into two uniaxial nematic phases. It was believed that it can happen mainly 
due to attractive van der Waals interactions, but there is another possibility of phase separation. 
It was shown that a system of two kinds of hard spheres can undergo a demixing transition 
[23] (for the exact formulation of the problem of miscibility for mixtures of spheres see also 
[24]). In the case of rod-disc mixtures the demixing is also driven by steric interactions. It 
was proven using Onsager theory that only at elongation of rods equal to 15 and discs equal 
to 1/15 stable biaxial nematic phase is possible in a narrow region in number density [25]. 
That inspired further research on the issues of phase separation. The Monte Carlo simulations 
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study of hard discs of elongation 1/15 and 1/20 and hard rods with elongations 15 and 20 [26] 
confirmed the findings from [25]. N B was found to be stable in a region of molecular elon- 
gations and fractions of discs and rods significantly limited by areas of separated, coexisting 
uniaxial nematic phases. The ultimate solution to the demixing problem was proposed in 2003 
[27]. In that study a so-called shape amphiphiles were proposed, built by covalently linking 
the rod-like molecule to the disc-like one. The resulting compound possessed the molecules 
of high average biaxiality, but so far even that approach failed to provide the biaxial nematic. 
A picture of the shape amphiphile is shown in Fig. 1.9. 



The binary mixtures of rod-like and disc-like molecules were believed to be good candi- 
dates for the discovery of a stable biaxial nematic phase. However, after the intensive studies 
it was found that most probably the system undergoes a demixing transition into two coex- 
isting uniaxial phases Nu + and N^_, before the transition to the biaxial nematic takes place. 
Even if Nb can be stabilized for certain molecular elongations and/or intermolecular potential 
strengths as suggested in [18, 20, 25, 26], the region of stability against demixing has to be 
very narrow. In practice this prevents the observation of the thermotropic biaxial nematic in 
binary mixtures. Also, on the experimental level it is extremely hard to realize the suggested 
model where the molecules of different types are more attracted to each other than to their 
own kind. It is so because the oblate nematogens are significantly different than their prolate 
counterparts, and separated uniaxial phases are favoured over the biaxial mixture. 

In the mixtures the interaction between rods and discs is the source of biaxiality. Alterna- 
tively it can be introduced by taking into account a system consisting of molecules of the same 
type with a degree of shape biaxiality and interacting via a biaxial potential. In the subsequent 




(a) Structure of shape am- 
phiphilic compound from [27]. 



(b) Schematic picture of the connection 
of a rod to a disc. 



Figure 1.9: The amphiphilic molecule. 
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sections we take a look at the approaches that follow this path. 

D 2 h - symmetric molecules 

Straley proved that the simplest deviation from the uniaxial molecular symmetry, i.e., in- 
clusion of D 2 h - symmetric molecules (see Fig. 1.5), leads to a stable biaxial nematic [10]. 
Furthermore, he showed that at a certain shape anisotropy a direct transition from isotropic 
liquid to can be observed at an isolated point. This so-called Landau point (also known 
as bicritical, or self-dual) marks the place where the three nematic phases N[/_, N[/+, N#, 
and isotropic state "meet". An example of such place on the phase diagram can be seen in 
Fig. 1.7. At that point the biaxial nematic is most stable, in the sense that there takes place 
a direct second order I so - Ng transition. It significantly changes the topology of a phase 
diagram, therefore it is of interest. In 1989, Mulder published a paper on isotropic sym- 
metry breaking bifurcations, where a condition for this point was derived [28]. Then it was 
tested for Straley model and hard spheroplatelet using Onsager approximation. Further studies 
on spheroplatelets and systems of biaxial ellipsoids confirmed the condition for the Landau 
point [29] where Density Functional Theory and the so-called smoothed (or weighted) den- 
sity approximation [30, 31] were used. In that approach the free energy part that includes the 
contribution from intermolecular forces is taken as the free energy of isotropic fluid calculated 
at "smoothed density", which is related to the actual density by a formula involving the pair 
excluded volume. Landau points were identified for a range of molecular dimensions, confirm- 
ing the earlier findings. For hard ellipsoids the Landau point was found to occur for molecules 
with axes fulfilling square root rule: a x = ^Ja y a z (for definitions of a x , a y , a z see Fig. 1.5) 
which was called a self-dual geometry. Furthermore, the transition density at this geometry 
was found to be significantly lower than for close-packing, which suggests that the biaxial 
phase realisation is possible. These findings were confirmed by Monte Carlo simulations of 
hard ellipsoids with full translational and orientational freedom [32]. The biaxial nematic was 
found to be most stable near the self-dual molecular geometry, in agreement with previous 
studies, also the existence of Landau point was confirmed, and I so, N[/+, Nj/_, and N B states 
were identified. In 1997, the simulations of D 2 h - symmetric hard ellipsoids presented in [32] 
were studied in more detail [33]. Again the position of Landau point was confirmed to occur at 
the self-dual geometry. The isotropic, two uniaxial nematics, and biaxial nematic were found. 
Additionally it was found that the first order I so - r% transition with increasing molecular 
shape biaxiality becomes weakened; even for small deviations from uniaxial symmetry the 
discontinuity associated with first order transitions is considerably reduced, in agreement with 
earlier indications [34]. 
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The above studies concerned mainly the hard core interactions, which are practical, yet 
are known only to accurately reproduce the quantitative properties of the phase diagram. A 
more complete study, including the attractive dispersion forces, on the existence of the biaxial 
nematic was presented in the year 2000 [35]. Berardi and Zannoni showed that Monte Carlo 
simulations of a system of biaxial molecules interacting via a biaxial version of the Gay-Berne 
potential (developed 5 years earlier [36]) give a stable biaxial nematic. The authors studied 
a single molecular geometry in the rod-like region, with a number of interaction strength pa- 
rameters sets. Surprisingly, Nb was found only when the lateral attractive forces for biaxial 
ellipsoids were dominant, that is, the configuration where the molecules are facing their sides 
was preferred. 

The research on biaxial molecules confirmed the mean field predictions of Freiser [5] and 
Straley [10]. Firstly, both Density Functional Theory and Monte Carlo simulations predicted 
that purely repulsive, steric forces between molecules of D 2 h symmetry give rise to the stable 
biaxial nematic phase. Furthermore, it was proven that for certain molecular dimensions there 
exists a Landau point, where system undergoes a second order transition from isotropic phase 
directly to the biaxial nematic. All the studies confirmed that this is an isolated self-dual point; 
so far there are no microscopic phase diagrams with a line of I so - Nb transitions, although 
Landau theory can predict such behaviour [37]. 

The theoretical approaches predicted the existence of stable biaxial nematic in the system 
of rigid molecules of D 2 h symmetry. There exists a class of mesogens behaviour of which due 
to their molecular structure cannot be accurately modelled by systems of such objects. Those 
materials, known as bent-core systems, exhibited a variety of new phases, and also proved to 
be of importance in matter of thermotropic biaxial nematic. In the following section we review 
some of the properties of these systems. 

Bent-core systems 

In the so-called bent-core systems, the molecular structure resembles that of a boomerang, 
or a banana. They behave like entities of C 2 symmetry (see Fig. 1.10(b)), and were found to 
exhibit the biaxial nematic phase. The first theoretical study aimed at understanding of the 
biaxial ordering in bent-core systems was based on a model of two, hard, rod-like sphero- 
cylinders joined at the ends at a given bend angle (also called opening angle). Onsager theory 
was applied together with bifurcation analysis, and the phase diagram in plane of reduced 
density and bend angle was presented. It showed the isotropic phase, two uniaxial, and biaxial 
nematic. At the angle of 107° the Landau point was found at the value of reduced density 
accessible in experiment [38]. Similar results followed from the mean field analysis of a lat- 
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(a) Diagram predicted by 
mean field analysis of a lattice 
model of V-shaped molecules 
as presented in [39]. 



(b) Model bent-core 
molecule with opening 
angle 7. Dashed line 
represents C2 symmetry 
axis. 



(c) Bifurcation diagram for 
bent-core molecules com- 
posed from two hard sphero- 
cylinders from [38]. 



Figure 1.10: Predictions of phase diagram for model bent-core molecules with bend 
angle 7. 



tice model composed of V-shaped molecules [39]; both diagrams are shown in Fig. 1.10. The 
issue of a possible biaxial nematic and phase sequence in bent-core systems was also widely 
addressed by simulations. In 1999, a Monte Carlo study of a system of bent-core molecules 
composed of two hard, oblate spherocylinders was presented [40]. The elongation of the 
resulting banana-like molecule for the bend angle of 180° was chosen to be 4:1. A slight 
deviation from the rod-like shape resulted in appearance of uniaxial nematic for the opening 
angle between 160° and 170° in higher densities followed by smectic-A phase. Smaller angles 
between arms destabilized the nematic phase and system crystallized to a biaxial solid. For 
angles of 120° and 110° molecules formed interlocking pairs and thus orientational ordering 
was not seen. Also the case of elongation 18:1 (in the rod-like limit) was studied for bend an- 
gle of 108°; the simulations in the biaxial nematic proved it to be mechanically stable against 
isotropic liquid. No other sign of the spatially uniform biaxial phase was found. The inclusion 
of Gay-Berne interacting parts as model arms of the bent-core molecules also failed to pro- 
duce the elusive phase. For the model of two Gay-Berne arms each elongated in proportion of 
3: 1 Monte Carlo simulations for bend angle of 140° showed that isotropic phase was followed 
by uniaxial nematic and then by smectic-A [41]. For the opening angle of 170° the isotropic 
liquid lost stability at the expense of smectic-A, and no uniaxial nematic was seen; it appeared 
for lower angles of 160° at a single point [42], where also tilted smectic-B phase was seen. No 
alignment of the steric dipole axis was found [42] . 

The issue of polar structures, which may follow from the steric and electric dipoles present 
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in the bent-core molecules, was addressed by molecular dynamic simulations for a two Gay- 
Berne molecule bent-core model [43] where the electric transverse dipole was introduced in 
the plane of the opening angle. The dipole-dipole interactions stabilized the smectic-A and B 
phases at the expense of the uniaxial nematic. The influence of the location of the dipole on 
phase sequence was also studied for a model of the banana constructed from three Gay-Berne 
molecules by Monte Carlo simulations [44]. When the dipoles were localized on the arms 
of the molecule, the uniaxial nematic was observed. For a central dipole along 0% symmetry 
axis, the isotropic phase lost stability to smectic structures with lowering of the temperature. 

None of the above simulation studies detected any sign of the biaxial nematic. Some proof 
of biaxial ordering was presented in the atomistic simulations study [45], however, the degree 
of order was low. Only a simple lattice Lebwohl-Lasher model [46, 47, 48] exhibited a stable 
biaxial nematic in Monte Carlo studies. It was found that in the system of V-shaped molecules, 
the increase in asymmetry of the arms shifts the Landau point to lower angles [46]. Similar 
behaviour was observed when molecular flexibility was introduced by allowing the bend angle 
to vary [47]. The dipole-dipole interactions were also studied, and were found to lead to a new 
topology of a phase diagram, namely the line of Landau points in the range of opening angles 
at a certain dipole strength was observed [48]. 

Despite the predictions of the mean field and Onsager approaches for the bent-core mole- 
cules [38, 39], the simulations [40, 41, 42, 44] were unable (with the above exceptions 
[45, 46, 47, 48]) to find a stable biaxial nematic phase. Biaxiality was exhibited only in 
crystalline structures, and in the smectic phases. Usually the latter dominated the phase dia- 
grams obtained in the simulations. Those approaches also indicated that the polar long-range 
nematic order due to the presence of the steric and/or electric dipoles is less stable than non- 
polar nematic and smectic structures. 



1.3.3 Unsuccessful experimental attempts until year 2003 

The predictions of Freiser and Straley of the early 70s and later theoretical and simulation 
studies suggested that the biaxial nematic can be stabilized in a system of molecules pos- 
sessing "sufficiently" biaxial shape. That inspired a series of experimental studies (following 
the discovery of Yu and Saupe in 1980 [11]) focused on the synthesis of compounds with 
molecules that possessed a structure which would bring enough biaxiality to the system to 
make the observation of the fhermotropic biaxial nematic phase possible. Fig. 1.11 shows the 
approximate molecular shapes that were advanced as good candidates for the realisation of 
N B phase [49]. All of those were tested experimentally and apart from one, which was found 
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(a) Biaxial ellipsoid (b) "Dumbbell" - like (c) "Pipette" pro- (d) "Lollipop" motif 
motif considered in shape studied in [54, posed in [50] from [53] 

[51,52]. 55] 

Figure 1.11: Schematic representations of molecular motifs tested in unsuccessful 
attempts of finding the biaxial phase [49] . 

to be uniaxial [50] (Fig. 1 . 1 1(c)), were claimed to form the fhermotropic biaxial nematic phase 
[51, 52, 53, 54, 55], mainly on the basis of the optical observations of textures. They were all 
reinvestigated by means of 2 H NMR and the nematic phases were found to be actually uniaxial 
or the phase exhibited a too small degree of biaxiality to consider it to be truly biaxial. 

The attempt to stabilise the thermotropic biaxial nematic based on the molecular design of 
low-mass compounds failed, but the search continued. Finally, after over thirty years of strug- 
gle, at the beginning of year 2004 came the long-expected discovery. Next section is devoted 
to a brief description of experiments showing stable biaxial nematic phase in thermotropic 
materials. 



1.3.4 The discoveries after 2003 



The first convincing report on the thermotropic biaxial nematic phase was presented by 
Severing and Saalwachter in March, 2004 (the paper arrived within the editors four moths 




Figure 1.12: Side-on attachment of mesogen to a polymer chain [56], giving rise to 
the biaxial nematic phase. 
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(a) (b) 



Figure 1.13: Structure of two bent-core compounds; (a) Bent-core molecules studied 
in [58], opening angle in biaxial nematic phase was estimated to be 140°. (b) Bent- 
core molecule with bend angle of 90° (after [61]). 



earlier) [56]. The authors studied a system consisting of molecules composed of rod-like 
mesogen attached to a polymer (see Fig. 1.12). That connection of the low-mass nematogenic 
compound to the side-chain polyatomic structure resulted in the biaxial nematic phase of the 
attached mesogens. The discovery was confirmed by 2 H NMR. 

A few weeks later, in April of 2004, two independent reports were published, first one by 
the Kent group (arrived within editors in July of 2003) [57], second in October [58] (being a 
reprint and a confirmation of earlier, preliminary study [59]), both concerned the experimen- 
tal investigation of the thermotropic biaxial nematic phase in systems of low-mass bent-core 
molecules of similar structure, see Fig. 1.13. As can be seen, each molecule has a strong elec- 
tric and steric dipole moment along C 2 axis, and two flexible "tails" connected to the ends of 
arms. The findings of the biaxial order in the x-ray diffraction [58] were confirmed by the 2 H 
NMR experiments [57], and the bend angle in N B phase was estimated to be about 140°. A 
year later another report of a stable biaxial nematic order for a different bent-core compound 
was presented [60]. Nb was observed using x-ray diffraction and optical studies following 
N[/ + phase, and when the temperature was lowered further, a sequence of three smectic phases 
(smectic-C, X and Y) was detected. Very recently the bent-core mesogens with the opening 
angle of 90° were synthesised and showed to produce both the biaxial and uniaxial nematic 
phases [61]. All these results for banana-shaped molecules were surprising; however, earlier 
it had been indicated that bent-core systems can exhibit the elusive phase [62]. 

The third realisation of the thermotropic biaxial nematic phase came from systems of 
molecules which in their structure link the properties of two bent-core mesogens; the so- 
called tetrapodes. These compounds consist of four mesogenic molecules connected to a rigid 
siloxane core through siloxane spacers. In result the system forms the molecules resembling 
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Figure 1.14: Structure of the tetrapode compound for which the biaxial nematic phase 
was found as shown in [63], and its schematic shape. 

on the average a fiat platelet (Fig. 1.14). In 2004, Merkel, Kocot, and co-workers reported 
that materials composed of such molecules form a stable fhermotropic biaxial nematic [63]. 
The findings were based on the infra-red absorbance measurements, which allowed for order 
parameters to be calculated. It was shown that the system undergoes a second order transi- 
tion from the uniaxial nematic to the biaxial nematic phase. A significant degree of biaxial 
ordering was determined, the results were confirmed by conoscopy and texture observations, 
and the temperature dependence of order parameters was found to be in agreement with mean 
field predictions [64]. The phase behaviour of similar systems was later studied by means of 
2 H NMR and N B was also found [65]. 

The discovery of the biaxial nematic phase in systems of banana-like mesogens came 
as a surprise mainly because these compounds of low molar mass where known to produce 
mostly smectic phases, the behaviour which had been recovered in some simulations. Neither 
they were a newcomers to the world of liquid crystal physics; the first bent-core mesogenic 
molecule was synthesised by Vorlander in 1929 [66] and the banana-like material were a sub- 
ject of intensive studies since then. Those revealed a great variety of polar and chiral smectic 
structures, many of them seen for the first time thanks to the bent-core mesogens [67]. Since 
one of the threats for the biaxial nematic ordering is for a system to stabilise the smectic phase 
before the biaxial nematic, the materials best known for producing smectics would be the last 
place to look for a stable thermotropic biaxial nematic. However, it proved to be otherwise, 
and now it is clear that the bent-core and tetrapode systems are the only low-mass compounds 
where the thermotropic biaxial nematic was found to be stable so far. Therefore they are of 
interest. 

As we have mentioned, the molecules in uniaxial nematic phase are rotating relatively fast 
about "the longer" axis, which essentially causes the constituent of the rod-like nematogen 
to behave on average as uniaxial objects. The above examples of tetrapodes, bent-core and 
polymeric systems prove that effective hindering of that rotation may lead to the stable ther- 
motropic biaxial nematic. On the other hand, more than thirty years of research resulted in 
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four classes of substances (including the lyotropic system) giving rise to N#. They seem very 
different from each other (see Fig. 1.15) and from earlier proposals (see Fig. 1.11). Thus it is 
clear that the biaxial nematic order is not easy to induce. However, there probably exist certain 
factors that make the elusive phase more likely to occur. It is not straightforward to deduce 
them just by looking at Fig. 1.15 but we can investigate the models of biaxial nematic phase 
and try to localize some of those factors. This is the main goal of the present thesis, as will be 
described in the following section. 

1.4 Purpose of the thesis 

Present work is devoted to the investigation of some of the factors that influence the stabil- 
ity of the biaxial nematic phase. As we have seen, there are many unresolved issues concerning 
stabilisation of this liquid crystal state. In the case of banana-like molecules the theory dic- 
tates that N B is most stable, in comparison to Nu and I so, for the bend angle near 107°, which 
differs significantly from the experimentally estimated value of 140°. Also the influence of 
electrical dipoles, present in these systems, on biaxial nematic ordering is unclear. For D 2 h 
- symmetric molecules the hard interactions provided the stable Nb, however, for the more 
complex but also more realistic biaxial Gay-Berne model only one Monte Carlo study was 
presented, and biaxial nematic was stable only for one set of potential parameters. We are 
going to address those and other issues of stabilisation of Ng phase in a more systematic way, 
namely by taking into account models of bent-core molecules and studying the biaxial Gay- 
Berne interaction. We aim at finding approximate phase diagram and identifying molecular 
factors responsible for appearance of the biaxial nematic order. Much effort is devoted to the 
investigation of Landau points 2 , since there the action of microscopic model parameters on 
macroscopic biaxial nematic ordering is strongest, and its stability is enhanced. These points 
are indicated as the best places to look for a stable Nb. 

Up to date there are only few theoretical treatments that allow to predict a phase diagram of 
liquid mesophases from first principles. One of them is the Local Density Functional Theory 
(DFT), which is known to give results of reasonable agreement with simulations and experi- 
ment for a variety of liquid crystalline models. Yet to our knowledge it has not been employed 
in a systematic way to the systems with biaxial nematic phase. It is of interest to see how that 
theory works in this case. In majority of cases studied, we have decided to use this approach 
in low-density limit (also called second order virial approximation), and analyse the equation 

2 In present work the Landau points coincide with self-dual points obtained from duality transformation of 
L = 2 model, so the both names refer to the same physical entity. 
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(a) Side-on attachment of (b) Model banana-like (c) Schematic representation 
mesogen to polymer. molecule. of tetrapode (side view). 

Figure 1.15: Schematic structure of molecular motifs of the compounds that gave 
rise to the thermotropic biaxial nematic. 

for equilibrium by means of bifurcations to obtain behaviour of a system close to transition 
point. Both methods will be described in the following chapter. It is a well known fact that this 
approach overestimates transition temperatures, therefore we will not aim at the exact values 
of critical parameters, but rather try to identify possible topologies of the phase diagrams 3 . 

Firstly, we consider an effective model with three coupling constants based on Straley's 
[10] generalization of the orientational part of the Maier-Saupe potential [6], and present the 
exact bifurcation diagram in mean field, including spatially uniform phases. It serves as an 
introduction of some general concepts concerning studies on nematic order, in particular the 
case when the coupling constants are interpreted as coefficients of expansion of the pair direct 
correlation function used later. The derived diagram shows Landau (self-dual) and tricritical 
points. Next, the mentioned earlier biaxial potential of Berardi, Fava, and Zannoni [36] is 
studied using the low-density approximation. This model is the only soft interaction involving 
translational degrees of freedom for which the biaxial nematic phase was discovered in sim- 
ulations [35]. Therefore it presents a possibility of comparing our results with those obtained 
from Monte Carlo and then investigating the potential parameter space further. We begin with 
calculation of the bifurcation diagram for parameters sets studied in the simulations, and also 
compare the method of minimisation of the Helmholtz free energy with bifurcation for uni- 
axial - biaxial nematic transition. Then we address the issue of the Landau point, and locate 
it by varying the shape and energy biaxiality. Finally, we study the bent-core models. In that 
case we pursue the shape related effects as well as the influence of dipole-dipole interactions, 
in search for the decisive factors that stabilise the biaxial ordering. We consider the models of 

3 Most of the diagrams will be presented at bifurcation point, therefore will be approximated, possible phase 
diagrams. 
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Gay-Berne interacting parts by choosing the arms of bent-core molecules to be soft, prolate 
ellipsoids. We study the cases of two and three uniaxial arms with elongation 5:1, and also 
take into account two biaxial ellipsoids as the model arms. Then, we turn our attention to 
the elastic issues and study the bulk elastic free energy density in absence of chiral order by 
calculating the set of elastic constants in the biaxial nematic phase for the biaxial Gay-Berne 
interaction. We show the temperature dependence of the constants in prolate molecular regime 
and at Landau point. Finally, we present the bifurcation equations with the inclusion of the 
transitions involving orthogonal smectic phases of uniaxial and biaxial symmetry. We also 
briefly address the issue of apparent importance of lateral interactions in the stabilization of 
Nb [35], by presenting temperature dependence of order parameters including biaxial and uni- 
axial smectic-A phases for two sets of potential parameters, calculated for biaxial Gay-Berne 
in L = 2 model. We show how the temperature range of biaxial nematic, limited by biaxial 
smectic-A, gets wider for the case of strongest lateral forces. 

This work is aimed at the investigation of the way transition point to the biaxial nematic 
phase is affected by shape of molecules and parameters of the intermolecular forces strengths 4 . 
These microscopic quantities are related to the macroscopic properties of the system by Local 
Density Functional Theory, and the transition temperatures and densities are approximated 
by means of bifurcation analysis of the self-consistent equation for equilibrium one-particle 
distribution function. In that way we obtain a possible phase diagram parametrized by the 
constants entering pair potential, which can be used as a guidelines of the role different mi- 
croscopic quantities play in the stabilisation of biaxial nematic. 

The thesis is organized as follows. Firstly we present the general introduction to the Lo- 
cal Density Functional Theory, which is particularly useful in the studies of phase transitions. 
Next we describe a general bifurcation scheme for finding the point where a new, lower sym- 
metric solution branches of the reference one of higher symmetry, together with methods of 
finding the Landau and tricritical points. In subsequent chapters the three models describing 
the biaxial ordering are studied. Firstly, the model of three coupling constants in a general 
expansion of D 2 h - symmetric pair interaction is investigated in mean field approximation. 
Then, we take into account the biaxial Gay-Berne potential. We start by comparing the DFT 
results with Monte Carlo simulations and with minimisation of free energy, and continue to 
investigate the interaction parameter space further; the roles of shape and energy biaxiality 
are studied while acting separately and simultaneously, and the behaviour of Landau point is 
considered. Next, we take into account the bent-core molecules models, and pursuit the issues 
of shape and dipole-dipole interaction. In concluding chapters we briefly address the issues of 

4 In the models used in the present work, both the molecular shape dimensions and forces strength parameters 
are all incorporated into pair interaction potential. 
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elastic constants for biaxial Gay-Berne model, and comment on uniaxial and biaxial smectic- 
A phases, including the calculation of the order parameters in L = 2 model, and the derivation 
of bifurcation equations for additional transitions. 




Chapter 2 
Local Density Functional Theory 



The Local Density Functional Theory is a tool that presents a way for determining the equi- 
librium properties of a system from first principles. Once applied to liquid crystal models, 
it allows to determine the phase diagram. In the present work, using this method, we obtain 
the densities and temperatures for phase transitions in the systems with stable isotropic, uni- 
axial and biaxial nematic phases. Current chapter is devoted to a systematic approach to the 
Local Density Functional Theory. Also a particular form of the theory used in this thesis is 
introduced. 

2.1 General formulation for one-component systems 

In this chapter we introduce the theory that is used as a tool in the analysis of models with 
the biaxial nematic phase. We consider the liquid crystal mesophases by postulating pair in- 
termolecular potential and seek the phase diagram in density-temperature plane. In order to do 
this, we need to establish a connection between the microscopic parameters associated with 
molecular dimensions and interaction strengths, present in a model, and macroscopic system 
properties. We are interested in the behaviour of the system at equilibrium in the vicinity of 
phase transition point. Local Density Functional Theory (DFT) provides a way of systematic 
introduction of intermolecular forces. It can be used to determine the conditions for equi- 
librium, and to establish the connection between microscopic and macroscopic parameters. 
Furthermore, it allows to derive exact formulas for the point where a new structure emerges 
from a given equilibrium state. 

In this section we present the main concepts of DFT and derivation of the most important 
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formulas, following [68, 69, 70] l . The results were published by us in [71]. Although the 
approach can be easily modified to include molecular flexibility or to model polyatomic sys- 
tems, presently we restrict ourselves to the case of rigid molecules, orientation of which in 
euclidean space can be described by supplying the right-handed, orthonormal tripod of vec- 
tors associated with a given molecule, and the position by the vector linking to the molecular 
centre of mass. 

Let's consider a system of rigid, identical, biaxial molecules in grand canonical ensemble, 
where volume V, temperature T, and chemical potential fi are fixed. State of each molecule is 
described with respect to the global laboratory frame by the position Fj of centre of mass, and 
orientation f2j of reference frame associated with the molecule with respect to the laboratory 
system, and parametrized, e.g., by three Euler angles a^(5 i: ji [72, 73]. In the following we 
will use the notation x t = (r^fij) and J dxi = J dri J Q 27r dai §_ x dcos{(3i) d^ to rep- 
resent the degrees of freedom of the molecules, and available positional-orientational phase 
space. The system considered can be described by the Hamiltonian of the following form: 

H N = T N + U N + V N , (2.1) 

where T N , Un,Vn are the kinetic energy, the potential energy, and the external field interac- 
tion, respectively. They read: 

N ^ 2 i N 
i=l i=l 

U N = U(x 1 , . . . ,x N ) , (2.2) 

N 

vn v ext {xi) , 

8=1 

where N stands for the number of molecules, p, for the momenta, m for mass, uj\^ U2,%, 
denote the angular velocities, I l7 I 2 ,l3 are the principal moments of inertia, V ext (xi) is the 
external potential, and U(xi, . . . , x N ) is the potential energy. As we can see we are not mak- 
ing any assumptions concerning U, especially we do not yet introduce two-body terms. The 
class of systems considered includes all systems where the kinetic energy can be integrated 
out. Therefore, we will neglect the momentum variables, since the corresponding contribution 
does not affect the phase diagram. 

Each microstate is realized with a certain probability; it is described by the distribution 
function f(x%, . . . , x^), which is the probability density of a configuration where iV molecules 

'it is probably worth noting that a certain level of generality is maintained here. We believe it will stress out 
the flexibility and tractability of the approach. 
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are in a given state, taken between [xi, . . . ,xn] and [x± + dx±, . . . ,xn + cIxn], therefore 
called N-particle distribution function. Following Evans [70] we consider a functional S7 
of / (for clarity we will drop the arguments), defined as 

n [/] = Tr [{H N -fiN + /T 1 In /) /] , (2.3) 

where fi^ 1 = k B T, with Boltzmann constant k B . We used the following notation for a "clas- 
sical trace": 



°° 1 f 

Tr (^) = h3N N \ / dXl ■ ■ ■ dx NM X l, 



where h is the Planck constant. Clearly, the functional fl [/] for an equilibrium distribution 
f eq reduces to the grand potential 

n = n[/ e? ]=^- 1 lnH, (2.4) 

where the grand partition function 

E = Tr{exp[-P(H N -(jiN)]} , (2.5) 

and 



f eq = E- 1 exp [-p (H N - (J.N)] . (2.6) 

As usual, f eq is normalized with respect to the classical trace: Tr(/ e? ) = 1, and the ensemble 
averages are calculated with the help of f eq , namely (A) = Tr(f eq A). 

Let's consider an arbitrary distribution A of unit trace (Tr(A) = 1). The functional Q [A] 
can be written as: 

Q [/i] = Tr { [-/T 1 (In S + In f eq ) + /T 1 In A] A } 

= ft [/ e? ] + /5Tr (A In /i - A In f eq ) (2J) 

= n [f eq ] + /3Tr In A^J . 

Since m(A<j/A) < f eq /f\ — 1, where equality holds for A = f eq , we can see that the last 
term is positive-definite. It represents the well-known relative entropy S re i of Kullback and 
Leibler, introduced in 1951 [74]. From definition S re i[f 1 /f 2 ] = — fc^Tr [A In (A//2)], and 
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always S re i[f\/ f?\ < 0. So we can conclude that for any normalized / ^ f eq we have 



n [/] > n [f eg ] . 



(2.8) 



It means that there exists a variational principle, which states that the global minimum of 
S7 [/] is attained when / coincides with equilibrium distribution f eq . Another key quantity is 
a functional T [/], which for / = f eq plays a role of the Helmholtz free energy of the system 
in the absence of external field. It is defined as follows: 



and is traditionally called the intrinsic Helmholtz free energy functional, while the total Helmholtz 
free energy reads F = Tr (f eq V N ) + J 7 [f eq \. 

The knowledge of the N-particle distribution gives a complete information about the be- 
haviour of a system. Fortunately, in the description of equilibrium properties we do not need to 
know all details of f eq . Actually, it turns out that the equilibrium behaviour can be determined 
by noticing that both the grand potential Q and the Helmholtz free energy T are functionals 
of only one-particle distribution g(x), since / is also a functional of g(x) [70]. Even more, it 
can be proven that both Q and T are unique, although unknown explicitly, functionals of g[x) 
[70]. Using this general theorem we can proceed to construct approximation schemes for Q 
and T, we will now describe the most common ones. 

We start by introducing the one-particle distribution function g(x), which in equilibrium 
is denoted by g eq (x) and defined with the help of the microscopic density operator 



T[f]=Tx[{T N + U N + (3- 1 \nf) f] 



(2.9) 



N 




(2.10) 



i=l 



as an equilibrium average of the above, over the grand canonical ensemble: 



Qeq(x) = (g (x)) = Tr[f eq g (x)] . 



(2.11) 



g eq (x) is normalized to the average number of particles (N), 




(2.12) 
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Now we rewrite (2.3) with the help of (2.9), and obtain [69, 70]: 



17 [g] = T [g] - J dx g(x) + J dx g(x)V ext (x) 



(2.13) 



which for equilibrium one-particle distribution function reduces to the grand canonical poten- 
tial: 

ftiSeq] =3 r lQeq} + J dxTx[f eq g (x)]V ext (x)-n J dxTx[f eq g(x)} 



=Tr 



J dxg (x) V ext + T N + U N + (3 1 In f eq - /x J dx g 



x 



f, 



eq 



(2.14) 



=Tr (H N -fiN + In f eq ) f eq = U. 



T [g] and 17 [g] are related by the Legendre transformation, since fl [g] can also be treated as 
a functional of ip — jx — V ext , 



T [g] — f2 [ip] + dxi/j (x) g(x) . 



(2.15) 



Since in the first place we are interested in equilibrium properties of a system, we need 
to determine the condition which can be used to find g eq (x). Keeping in mind that 17 [g] is a 
unique functional of one-particle distribution function, the variational principle (2.8) leads to 
the following (necessary) condition: 



5Q [g] 



5g (x) 



0, 



(2.16) 



e=0eq 



which according to (2.13) can also be written as 



8g [x\ 



T[g] 



dx 2 [V ext (x 2 ) - n] g{x 2 ) 



0. 



(2.17) 



Eq. (2.16) is a consequence of the Hohenberg-Kohn-Mermin theorem [68, 69] stating that the 
minimum of the functional 17 [g] is attained when g coincides with equilibrium one-particle 
distribution function. In other words, g eq can be found by minimizing (2.13), for which the 
equation (2.16) is the necessary condition. As we can see, the minimum of 17 [g] can be 
found as a minimum of T [g] with the condition of normalization of g(x), where the chemical 
potential plays the role of the Lagrange multiplier. So we turn our attention to the functional 
representing the intrinsic part of the Helmholtz free energy. We can immediately divide it 
into two parts: jF id [g] describing the ideal gas, and the excess part, T ex [g] that includes the 
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contribution from intermolecular forces; 

F [q] = Fid [g] + Tex [q] ■ (2.18) 

The ideal part is given by 

T id [q] = /T 1 J dxg(x) {In [Ag(x)\ - 1} , (2.19) 

where A = ^h l2 (3y{2 

7r) 6 m 3 /i/ 2 ^3 is the constant resulting from the integration over mo- 
menta. Using the above we can employ the variational principle expressed by Eq. (2.17) to 
finally get 

In {Ag(x)} = +P{H- V ext (x)} . (2.20) 

OQ \X) 

Since the yet unknown T ex [g] depends on the one-particle distribution function, the above 
formula becomes the non-linear, self-consistent equation for g(x). It reads 



g(x) = A^e^exp 



5 -^ + V ext (x) 

OQ (X) 



(2.21) 



For given 5T ex [g] /Sg(x), the above equation can be solved for g(x), e.g., in an iterative man- 
ner. The solutions coincide with minima, maxima, or saddle points of the grand canonical 
potential functional; the equilibrium distribution is the one that minimises Q [g] . For V ext = 
and for U n invariant under a global translation and rotation, Eq. (2.21) always possesses a triv- 
ial isotropic solution g{x) = const. The fact that b~T ex [g] /Sg(x) depends on g in a complex 
way makes the existence of many non-trivial, non-isotropic solutions of Eq. (2.21) possible. 

In order to seek the non-trivial solutions of the equation (2.21), we need to turn our atten- 
tion to the excess part of the Helmholtz free energy, T ex [g] . It is of course impossible to carry 
out the calculations without some approximations. A standard method to approximate T ex [g] 
is to perform a functional Taylor expansion about some reference state 2 described by g re f{x). 



2 Here we do not yet make any specific assumptions as to the form of the reference state or its relation to the 
other states that may became stationary for different temperature and/or density. 
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The expansion reads 



Fex [o] =Fex [Qref] + J dx 



8F e x [o] 



8g (xi] 



1 f , , 5 2 T ex [g] 



{q(Xi) ~ Qref(Xl)} + 
Q=Qref 

Q=Qref 



(2.22) 



The above expansion should be treated with caution, not only because its existence can be 
questioned, but also because it requires a detailed knowledge of the reference state, determina- 
tion of which in general case can be a great challenge in itself, even with some approximation 
for the derivatives of the Helmholtz free energy. Here we restrict ourselves to T ex which is 
analytical in g(x), and so the Taylor expansion (2.22) is assumed to exists and converge. This 
approach apparently works quite well in determining the phase diagrams, as can be verified 
by making a comparison with computer simulations (see, e.g., [75]). 

The idea of treating the free energy as an expansion in g(x) is similar in essence to the 
Landau-de Gennes approach [8]. Actually Eq. (2.22) gives the thermodynamical foundation 
for this theory. Also from the above expansion one can derive exact formulas for, e.g., the 
Oseen-Zocher-Frank elastic constants, assuming g re f{x) describes an undeformed state, as 
well as get an accurate description of the freezing transition (Ramakrishnan-Yussouff theory 
[76]). 

The expansion (2.22) is traditionally rewritten with the help of a set of the so-called di- 
rect correlation functions c„, for which T ex is by definition a generating functional. More 
specifically, they are introduced as functional derivatives of T ex [g] with respect to g{x): 



Ci(xi, [Qref]) = - P 



C 2 (Xi,X 2 , [Qref]) = ~ P 



SFex [g] 



Sg(xi) 
8 2 T ex [g] 



e=8ref 



8g(xi) 5g(x 2 ) 



8=Sref 



(2.23) 



t-n • • • j x n , [g\) 



5c 



n-l 



Sg(x n ) 







5 n T ex [g] 



Q=Qref 



5g (xi) . . .5g (x r 



8=Qref 



Interestingly, each c n is related through integral equations to the set of ordinary distribution 
functions {g(xi), £2(^1, x 2 ), ■ ■ ■ , g n { x \i ■ ■ ■ > x n)}- F° r example, the first function of the set 
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(2.23) can be connected to g(xi) as (see Appendix A): 

Cl (x u [q]) = {h- V ext (xi)} - /T 1 In {Kq{x x )} , (2.24) 

and c 2 can be related to the pair correlation function h 2 (x\, x 2 ) = g 2 (xi, x 2 )/ g(xi)g(x 2 ) — 1, 
by the Ornstein-Zernike relation: 

h 2 (xi,x 2 ) = c 2 (x 1 ,x 2 , [g]) + J dx 3 c 2 (x 1 ,x 3 ,[g])g(x 3 )h 2 (x 3 ,x 2 ) . (2.25) 

Now (2.22) can be written as 3 : 

00 „ 

Tex [g] = Tex [g re f} - ^ / tfai • • ■ dx n c n {xx, x n )5g(x 1 ) . . . 5g(x n ) = 

n=l J 

f-u+r'/^w^w- <2 ' 26) 

i y dx 1 dx 2 Sg(x 1 )Sg(x 2 )c 2 (x 1 ,x 2 ,[g ref ]) - . . . , 

where 5g(xi) = g(xi) — g re f(xi). With the definitions (2.23) and for V ext = 0, the self- 
consistent equation (2.20) can now be rewritten in a compact form as 

e( x ) = Zg 1 ex P { [q])} , (2.27) 

where Z g = J dxexp {ci(x, [g])} / (N) assures the normalization (2.12). We can conclude 
that the equilibrium one-particle distribution function is determined by the effective one- 
particle "potential": — /3 _1 ci(x, [g]), in a way specified by Eq. (2.21). In what follows we 
set V ext = 0, for this work is devoted to the liquids in the absence of external fields. 



2.2 Bifurcation analysis. Exact results. 

Information about the behaviour of a system close to a point of phase transition is stored 
in the self-consistent equation (2.27). Due to the structure of this equation, we can locate the 
points where a new solution branches off the reference one. This, by construction, takes place 
at a critical point or close to a real phase transition if the transition is of first order. We can also 
find sectors in parameter space where character of phase transition changes from continuous 
to first order - the so-called tricritical points [77]. In current section, we present a derivation 



3 We have assumed the equality of the chemical potentials of the reference and actual states. 
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of bifurcation equations, which are used to obtain the phase behaviour of biaxial nematics in 
subsequent part of this study. We also show the condition for bifurcation point to be a tricriti- 
cal point. 

As mentioned before, the solutions of the self-consistent equation (2.27) describe the lo- 
cal extrema and saddle points of the intrinsic Helmholtz free energy T [g] . The stable state 
is associated with the solution in the class of local minima corresponding to the global min- 
imum of T [g] . Remaining local minima are usually associated with metastable states, while 
the maxima are not realized as macroscopic states of a system. At sufficiently low density 
(or high temperature) T [g] has only one minimum corresponding to the isotropic distribution 
g(x) = const, which describes the unordered, isotropic phase. We expect that when the den- 
sity is increased (or the temperature lowered) the system undergoes a phase transition, which 
means that a new solution of Eq. (2.27) and a new local minimum of T [g] emerges. In cases 
studied here the symmetry group of high temperature phase contains all elements of the low 
temperature state. Therefore we say that the former has higher and the latter lower symme- 
try, and the relation between them is of group-subgroup type. The point where the higher 
symmetric solution loses the property of being the minimum of the free energy and becomes 
unstable with respect to the lower symmetric one which branches off is called a bifurcation 
point. For the first order transitions, it corresponds to a spinodal point, thus the bifurcating 
solution describes a metastable state. In practice, for the first order Iso - r% transition, it 
means that at the bifurcation density the minimum associated with the isotropic state changes 
to a local maximum of the free energy, and the only local minima left correspond to the uniax- 
ial nematic phases. For continuous phase transitions, the bifurcation point corresponds to the 
critical point, and the lower symmetric solution becomes immediately stable. Since present 
study deals with the weakly first order and second order transitions, the bifurcation analysis 
provides good estimates for the transition points. 

The choice of the class of phase transitions where the states are connected by group - 
subgroup relation and analysis of the solutions of the self-consistent equation that branch off 
the reference (higher symmetric) one constitute the so-called symmetry breaking bifurcation 
analysis. For the first time it was applied to liquid crystals by Kayser and Raveche for Onsager 
model of hard, long rods [78], and later was generalized by Mulder [28]. It is a numerically 
tractable method that allows to acquire some insight into the behaviour of the system close to 
phase transition. For the transitions studied in the present work it is a good choice, yet in gen- 
eral case it does not exhaust all the possible scenarios for the phase sequence. It does not tell 
which of the bifurcating states will be preferred. However, we can determine the symmetry of 
bifurcating state and character of the bifurcation. Present section is devoted to the description 
of this approach. 
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We begin with finding the exact equations for points where a given solution branches off 
the reference one, then continue to present the exact bifurcation equations for the spatially 
uniform phases and their specific form for Iso - Njj, Iso - Ng, and Nry - N# transitions. 
Next, we make a remark on a possibility of Landau points, and in the concluding section 
present formulas for the tricritical point. We follow the scheme presented by Mulder for the 
isotropic reference state [28] and later extended by us to the case of reference state of arbitrary 
symmetry [71]. 



2.2.1 General bifurcation equations 

Let's start by noting that a straightforward consequence of the expansion (2.26) is that the 
self-consistent equation (2.27) can be written as: 

6P(x) = Z re/ P re/ (x)Z; 1 exp jf> n A' n+ i (x, [*P]) J - P ref (x) , (2.28) 

where we have introduced average density p = (N)/V (with V standing for volume); 

and where 8P(x) = P s (x) — P re f(x) denotes a small deviation of the bifurcating P s (x) from 
the reference state P re f(x), and where 



K n+1 (x, 



If n 
[5P]) = - / c n+1 (x, x lt ...,x nt [8P]) J] SPixJdXi , (2.30) 

' J i=i 



with the normalization constant 



Z s = ^f- j dx P ref (x) eX P j^^n+X [P]) 

The normalization of P s (x) and P re f(x) implies that J dx5P(x) = 0. By construction, 
5P(x) = satisfies the equation (2.28) giving the self-consistent equation for P re f equiv- 
alent to Eq. (2.27). It is natural to seek the solutions 6P(x) ^ that branch off from the trivial 
one SP(x) = since they represent possible equilibrium states accessible from P re f state 
through a phase transition. As described in [28] and [71], these can be found by performing 
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the following expansions in powers of some arbitrary control parameter e: 

5P(x) = epx(x) + e 2 p 2 (x) + . . . , 

(2.31) 

p = po + e pi + e p 2 + • • • • 

Now inserting the above into Eq. (2.28), and comparing the terms of equal order in e we obtain 
the formulas fulfilled by p n at the bifurcation point. They read 



Pl(x) = Po Prefix) < K 2 (X, [pi]) - K 2 (x, [pi]) \ , 



p 2 {x) = po Prefix) \ K 2 (x, [p 2 ]) - K 2 (x, [p 2 ]) 



Pl P ref (x) < K 2 (X, [Pi]) ~ K 2 (X, [pi]) 



+ p 2 Q P ref (x) K 2 (x, \px]) K 2 (x, [pi]) - K 2 (x, [pj) 



(2.32) 



+ - \ K 2 {x, [ Vl ]f - K 2 (x, [px]) 2 \ + if 3 (x, [pi]) - K 3 (x, [pi]) 



p 3 (x) 



where all the averages in the above expressions are calculated with respect to the reference 
state: A = J dxA(x)P re f(x). Also, due to the normalization of one-particle distribution 
function, we have J dtip n (Cl) = and we can impose l/V J dxp n (x)pi(x) = 5„ 5 i 4 . Equa- 
tions (2.32) can be simplified further by recalling that the symmetry of the bifurcating state 
is lower than that of P re f and that they remain in the group-subgroup relation, that implies 
K n (x, [pi]) = and f dx P re f(x)pi(x) = 0. After substituting to the above, we obtain a set 
of hierarchical equations for pi (x) ,p 2 (x),...\ 

Pi(x) =p P re f(x)K 2 (x, [pi]) , (2.33) 

P2<» =P0 Pref(x)K 2 (x, [p 2 ]) + pi P re f(x)K 2 (x, [pi]) 

+ l -plPref{x) [K 2 (x, [ Pl ]) 2 - K 2 (x, [ Pl ]) 2 + 2K 3 (x, [ Pl ])] , (2.34) 



The first of the above is commonly known as the bifurcation equation, describing the branch- 
ing point to a new solution pi(x). Thus in the limit of small e, the low-symmetry state is 
P s (x) = P re f{x) + epi(x). The remaining equations determine the character of the bifur- 

4 It follows from the fact that Eqs. (2.32) are also fulfilled by p' n = p n + a n p\ (it expresses the free- 
dom of monotonic reparametrization of e in Eqs. (2.31)), where a n are arbitrary parameters. By setting 

l/V J dxp n (x)pi(x) = 5 n s we fix a n = 0. 
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cation (first order vs continuous), and also allow to obtain a condition where this character 
changes (tricritical or, generally, multicritical point). Interestingly, as we can see, once we 
know P re f (either we have postulated it or found it by solving the self-consistent equation), 
in order to determine the bifurcation point, it is enough to know the second direct correlation 
function. Higher order direct correlations contribute only to the subsequent equations, also in 
the hierarchical manner (in the first order equation only c 2 is present, in second order addi- 
tionally c 3 appears, and so forth). 

Possible scenarios for bifurcation are schematically depicted in Fig. 2.1. 

In the next section we show examples of applications of the formulas (2.28)-(2.34) for 
spatially uniform phases, and then continue to derive the bifurcation equations for the systems 
with isotropic, uniaxial nematic, and biaxial nematic phases. 




P = Po +e Pi +£2 P2. P = P„ + E 2 P 2 + E"P4 . P=Po + EPi + E 2 P 2 + - ■ 

p!<0, P 2 >0. P 2 <0, P 4 >0. p>0, P 2 >0, P 4 >0. 



(a) First order phase transition. (b) First order phase transition (c) Continuous phase transi- 

in vicinity of a tricritical point. tion. 

Figure 2.1: Generic bifurcation diagrams. Schematic representation of behaviour 
of leading order parameter near bifurcation point p . Path leading from (b) to (c) 
describes a change of character of phase transition as observed at tricritical point. In 
the reference state P re f the order parameter vanishes. 



2.2.2 Bifurcation equations in case of spatially uniform states 

The equations presented above are exact as long as we do not make any approximations 
for the direct correlation functions and the reference state. We can still make some general 
assumptions as to the form of the pair direct correlation and one-particle distribution function 
and acquire the form of the bifurcation equations that is easier to use. In case of the uniaxial ne- 
matic reference state, in c 2 we disregard the terms that depend on the director orientation. This 
assumption is consistent with low-density approximation used later, and is also equivalent of 
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c 2 describing poly domain sample where the director dependence is averaged out c 2 (xi, x 2 ) = 
j- J d 2 hc 2 (xi, x 2 , n). Also we assume that c 2 is invariant under the global rotation and the 
particle interchange operation, i.e., c 2 (a;i,:r 2 ) = c 2 (x 2 ,x 1 ), which imply that it depends on 
the relative orientations and positions of molecules. In present work we consider spatially 
uniform structures, namely isotropic, uniaxial nematic, and biaxial nematic phases. Thus, we 
take into account the bifurcations from Iso to Nu, from Iso to Ng, and from Nu to N#. For 
those transitions one-particle distribution function depends only on orientational degrees of 
freedom, i.e., P{x) = P(h), with / dfl P(Cl) = 1, P re f{x) = P re f{&), also^(x) = p;(fi), 
and pair direct correlation function c 2 (xi, x 2 ) = c 2 (f2i, f2 2 , ?i, r 2 ) = c 2 (Cl 1 f2 2 ,?i — r 2 ), 
where Q 1 Q 2 stands for relative orientation. With those assumptions we can present more 
specific and simpler form of Eq. (2.33). 

Close to bifurcation we can assume that pi(fl) can be expressed as 

P!(n)= J>g n Ag n (fi), (2.35) 

l,m,n 

where am,n are constant coefficients, and Am >n (fl) stand for real, orthogonal, symmetrized, 
linear combinations of Wigner matrices D$ >n (£l) [73, 79]. As described in Appendix C, 
A® n (f2) can be defined for arbitrary symmetry of phase and molecule, those correspond to 
the first and second of lower indices, respectively. In the above equation the summation goes 
over the physical and relevant for a given model ranges of indices; always —l<m,n<l, 
and, e.g., in case of D 2 j t symmetry all indices are positive and even. The coefficients am,n are 
unknown parameters which can be determined by means of bifurcation equations. Presently 
we proceed to find the bifurcation point p . The symmetry adapted functions AS,n(0) form 
an orthogonal base in the space of real functions, so we have the advantage of following 
orthogonality relations (for derivation see Sec. C.2): 

/Q 2 
dQ A« n (fl) A£? |B ,(«) = K A 6 ltV 5 m , m , 6 n , n < , Ni A = , (2.36) 

dn'A®, (n'-'n) ^(A') =n aa n^' aS iB (o) ^ s amyn , s ltl , , (2.37) 

o-,o-'e{-i,i} 

where N A ,n = (v / 2/2) 2+<5m0+5n '°. Using the above we can construct the equations for a,m,n 
fulfilled at bifurcation point; by inserting (2.35) to (2.33) and casting it on Am,n(0) we arrive 
at 

agn = {KaV 1 PoY, ai pi I ^^ , ^ e /(0)c 2 (n / - 1 n)AW(n , )AS )Tl (n) , (2.38) 

k,p,q 



36 



Local Density Functional Theory 



/-l", .o^ ,"'-1 



where c 2 (f2 ft) = J c/ 3 rc 2 (f2 f2, r). 

If we fix / = / in Eq. (2.35), i.e., if we choose one subspace of infinitely numerous set of 
subspaces labelled by the squared angular momentum index I, we can cast the above equation 
as 

4a} = { N L) ~* Po £ , (2-39) 

s 

where capital letters stand for two indices; {A} = m, n, and the summation over A means 
the summation over m and n in the manner described above, and where the elements of a 
bifurcation matrix £>® read 

^{Iub} = j ^^ / ^r e /(^)c 2 (^ 1 ^)Ag } (^)A« } (fi). (2.40) 

Eq. (2.39) is the eigenequation with eigenvalue (N AA /p ). The parameter p depends only 
on the elements of matrix u>^' and can be calculated as one of the roots of characteristic 
polynomial, namely from the solutions of the following equation: 



det 



0, (2.41) 



where 1 is the unit matrix. 

We need to decide which of the eigenvalues following from the above equation correspond 
to the bifurcation point. If we associate p with the density, then in the absence of coupling, 
when p = 0, we can only expect the isotropic phase. With the increase of density an assumed 
phase transition takes place, and in higher densities the bifurcation occurs. Therefore for non- 
zero po the relevant point where the reference solution loses the property of being stable is 
described by minimal value of po, thus the bifurcation point is chosen from the set of solutions 
of the Eq. (2.41) as the one corresponding to the lowest po- 

In general case, in Eq. (2.40) terms with higher angular momentum index than I can be 
present, but only due to the presence of P re f(fl). As can be seen from Eq. (2.37), the inte- 
gral over fl in Eq. (2.40) will leave only terms proportional to A^\(fl) (which is a conse- 
quence of the global rotational invariance of c 2 (ri ft)). It means that the leading contri- 
bution to the bifurcation from the pair direct correlation function will come from the part of 
c 2 (0 ft) ~ A^(f2 f2), which justifies the use of the expansion of c 2 in base of A^\(£l) 
functions. Only the presence of terms with A^ fc ](f2) with non-zero k in P re f(Cl) can generate 
in effect terms with higher angular momentum index in or - 1 (see, e.g. Eq. (C.19)). In the 
special case when P re f(fl) = const, the representations numbered by / decouple and uj^ 
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contains only terms of order I. It is studied in the following section. 



2.2.3 Results for D 2 h - symmetric model 

Presently we evaluate the above equations further and obtain the specific form of Eqs. (2.39)- 
(2.40) for the case of D 2h - symmetric molecules and phase. It means that the symmetry 
adapted base functions AS, n (Q) are defined for even and positive indices, and the summa- 
tions go over the relevant values of / > and < m, n < I, they read 

cr,cr'e{-l,l} 

where N™' n = (x/2/2) + m,0+ n '°, and where A$(fl) = 1. We also assume that for nematics 
we can truncate the expansion (2.35) at I = 2, i.e., that the leading contribution to bifurcation 
for spatially uniform phases I so, Nu, and Ng comes from coefficients am] n . It is consistent 
with the symmetry of these states and the choice of leading order parameters. 

We start with the special case when the reference state is isotropic which means that 

P ref (Cl) = i . (2.43) 

Then, the bifurcation matrix u)V> (2.40) contains only terms with given I, and the represen- 
tations labelled by Zs bifurcate independently (it is straightforward to check that this only 
happens when P re f(Cl) corresponds to isotropic state). 

Since, as we have mentioned above, the integral J dfl c 2 (^l Q)A m \ n (Q ) in (2.40) leaves 
in u>® only the terms with I, in order to find p , we can postulate the expansion of c 2 (f2 O) 
in the symmetry adapted base: 

c^n) = Y ^AgJ^O) , (2.44) 

where 

d(^ 1 ^)c 2 (0" 1 0)A«„(^ 1 ^), 

(2.45) 



c (0 = 2 J±1 / 
8tt 2 / 
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and where due to the particle interchange symmetry we have c$, n = c„, m . Now, using orthogo- 
nality of the base functions (2.36)-(2.37), and inserting the expansion of pair direct correlation 
function (2.44), Eq. (2.39) can be expressed in a simple form with I = 2, namely 



oL'L — ^ ] C q,n a rn,q • (2.46) 

96{0,2} 



iii.ii 



Above equation is a 4x4 eigenproblem. We recall that the first of lower indices of the base 
functions corresponds to the symmetry of phase. In this sense {a 2 Q, % 2 } are associated with 

(2) (2) 

uniaxial and {a 2 q, a 2 2 } with biaxial state. Another consequence of (2.37) is the fact that those 
two sets bifurcate independently, and in effect we get a 2x2 eigenproblem. It is so because 
c 2 (f2 f2) expressed by (2.44) gives the splitting of the space into subspaces labelled by 
angular momentum index / and phase related index m through the following term in (2.40): 
J dfl c 2 (f2 f2)Am,n(r2 ) ~ J2k^mk(^)- So, for isotropic reference state, not only the 
representations numbered by angular momentum index bifurcate independently, but also the 

phases of different symmetry branch off in separate subspaces. Furthermore, it is clear that 

(2) 

Eq. (2.46) is the same for uniaxial and biaxial a m \n- It can be easily solved, and we arrive at 
the following expression for bifurcation point: 

po = ■ 10 == • (2.47) 

,(2) , (2) _ L / (2)\ . / (2) J2) N 



C 0,0 + C 2,2 y^y C 0,2J + ^ C 0,0 C 2,2 

This equation describes the location of bifurcation point for transitions where the high- 
symmetry phase is the isotropic liquid and lower- symmetric state possesses spontaneously 
broken orientational symmetry 0(3) to uniaxial or biaxial D 2 h symmetry group. As we 
can see, the bifurcation point is fully determined by coefficients of the expansion of pair direct 
correlation function in the subspace of angular momentum index I = 2, and due to Eq. (2.37) 
depends on all of the coefficients Cm, n . Above expression is equivalent to Eq. (16) from [10]. 

Situation changes qualitatively when the reference state is taken to be a stable uniaxial 
nematic phase. In this case we need equilibrium, uniaxial P re f, to determine it we solve the 
self-consistent equation (2.27). The one-particle distribution function is taken to be of the 
following form 

P (") = E <A8 B >A£ B (n) . (2.48) 

l,m,n 

where the coefficients (A™ >n ) are calculated as the averages over the ensemble, and identi- 
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lied with order parameters (for the explanation of this model see Appendix D). Using the 
orthogonality relations (2.36) we can find that 

(A£„> = J dAP(A)Ag B (A). (2.49) 

From the above we choose the set with / = 2 as the leading order parameters, which is 
consistent with the symmetry of states in consideration. As we have mentioned, the first of 
the lower indices in Am in (Q) is identified with the symmetry of a phase, while the remaining 
one with a molecule. Along this convention (Aq 2 q) and (Aq 2 2 ) are associated with uniaxial 
nematic, and (A 2 2 q) and (A 2 2 2 ) with biaxial nematic. They all vanish in the isotropic phase. 
With apropriate choice of the moleuclar and lablatory axes, in r% ( Aq 2 q) ^ while (A^q) = 
(A22) = 0, for N B (Aq 2 q) 7^ 0, (A22) 7^ 0. The interpretation is clearer once we explicitly 
write the set of symmetry adapted functions Am )Tl (f2) for I = 2; the uniaxial functions read 



A$(n) =P 2 (cos{(3)) = \ + ~ cos (2(3) = - 1 - + ~ (b 3 -1 ' 



13 



(2.50) 



= ^Y cos (2 7) sin (/3 ) 2 = ^{(b 1 .i 3 ) 2 - (b 2 • 1 3 ) 2 } , 
where P2(x) = | (3 x 2 — 1) is the second Legendre polynomial, and the biaxial functions are 

Ag(O) = ^-cos (2 a) sin {(3) 2 = ^ | (b 3 ■ li) ' - (b 3 • 1 2 ) ' } , 

Ag> (fi) = - [3 + cos (2 (3)\ cos (2 a) cos (2 7) - cos (/3) sin (2 a) sin (2 7) = ( 2 - 51 ) 



'+(b 2 -i 2 ) 2 -^(b 3 -i 3 ) 2 -i 




Figure 2.2: Two orthogonal, right-handed tripods of unit vectors corresponding 
to the frames associated with laboratory ({li,l 2 ,l 3 }) and molecular (body) frame 
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where -Tlx , 1 2 , 1 3 J> and |b 1 ,b 2 ,b 3 | stand for right-handed, orthogonal tripods associated 
with laboratory (director) and molecular (body) frame, respectively (see Fig. 2.2). Now it 
is clear that (Aq 2 q) measures the degree of order of chosen molecular axis with respect to the 
laboratory-fixed direction, while ( A 2 2 ) reflects the ordering of complete tripod of vectors of 
molecular and laboratory frames. It is straightforward to check that — | < (Aq 2 q) < 1, and 
— 1 < (A 2 2 2 ) < 1. We can also see that when we choose 1 3 as the director, then the oblate ne- 
matic order N;y_ with respect to molecular axis b 3 entails in ideal case ( A q) = — | and ( A 2 ) 

~ A (2) (2") 

can be non-zero, while in ideal ordering of b 3 along 1 3 we have ( A q) = 1, (A 2 ) = 0. In 
these cases (A 2 5 1 ) = (A^ 2 2 ) = 0. On the other hand when \i is parallel to bi, i 2 to b 2 , and 1 3 to 
b 3 , i.e., we have ideally aligned N B , then (A 2 2 2 ) = (Aq 2 q) = 1 and remaining averages vanish. 

(2) A 

We can construct linear combinations of Am, n (f2) that satisfy the uniaxial symmetry restric- 
tions with respect to other laboratory axes, e.g., states a( Aq 2 q + y^A^) + 6(Aq 2 2 + \/3A^) 
and a(Ag 2 o — V3A^l) +6(Aq 2 2 — v^A^) describe phases uniaxial about i 2 and l x , respec- 
tively. (Aq 2 q) and (A 2 2 2 ) are referred to as leading, dominant order parameters. The remaining 

(2) 

two, as we can see from (2.50)-(2.51), (Aq 2 ) is sensitive to the orientation of the plane de- 
termined by two vectors associated with the molecule with respect to the laboratory chosen 
direction, while (A^) depends on the orientation of the plane determined by two vectors in 
the laboratory frame in comparison to the chosen molecular axis. 

In order to further work out the bifurcation equation (2.33) for the N[/ - N B transition, we 
need a stable reference phase. Thinking in terms of the expansion (2.48), we need to know 
the order parameters (Aq m ) for even / and for given temperature and density. Since now the 
reference state has non-trivial structure, we cannot expect that the bifurcation point obtained 
from (2.39) will depend only on order parameters with lowest possible angular momentum 
index 1 = 2. 

The reference phase can be determined with the help of the self-consistent equation (2.27). 
The equation can now be cast in the following form: 



P(Cl) = Z- 1 exp 



(2.52) 



where Z = JdQexp p J dQ c 2 (fi fJ)P(fi 



ensures the normalization. Once we insert 



the expansion of c 2 (f2 f2) (2.44) up to angular momentum index / = 2 to the above, we can 
write expl 
they read 



write explicitly the equations for order parameters using the definition of (Am,n) Eq. (2.49), 



(2.53) 



Local Density Functional Theory 



41 



where 



c£i(n)= Yl cg! B (A®>Ag(n) , 

m,ng{0,2} 



(2.54) 



m,ne{0,2} 



Setting (A^q) = (A^) = 0, i.e., c^J (f2) = 0, we can solve the self-consistent equations for 



,(2)\ 



.(2) 



(2) 

order parameters (2.53) interactively, provided we know the expansion coefficients c m ,n, and 
then use the order parameters to build the uniaxial reference one-particle distribution function 



Pre/fa) 



Z re )exp 



(2.55) 



where Z re f = J dfl exp pc 2u {tl) 

Having determined the reference state, using Eqs. (2.53)-(2.55), we can continue to find 
the equation for bifurcation point p . The eigenproblem (2.39) now involves a 2x2 matrix 



to 



(2) 

{a,b},{p,q} 



(2.56) 



where {a, b} e {2, 0}, {2, 2}, and {p, q} G {2, 0}, {2, 2}. Using the above with the expan- 
sion of c 2 (f2 ft) (2.44), orthogonality of base functions (2.36)-(2.37), and substituting the 
reference phase (2.55), we can cast the eigenequation (2.39) as 




A „(2) . r> „(2) o „(2) . 



Po 



(2) \ 

0,2 


M 


< <%l 


(2) , 
0,2 1 







(2.57) 



where 



A 2 = -A 



(2) 
0,2 



14 



A$ = 5o 



^2 = - + ^A 00 



70 



-A 



(4) 
0,0 



A (4) 



2V35 



t - t a( 2 ) + _a( 4 ) 

-^0,0 ^ or ^0,0 > 



5 7 u ' u 35 

and where the averages are calculated in the reference, uniaxial nematic state. The bifurcation 



42 



Local Density Functional Theory 



point can be calculated from the above eigenequation using (2.41), it reads 



2 




(2.58) 

Equations (2.47) and (2.58) describe the bifurcation scenarios for the system with three 
spatially uniform phases: isotropic, uniaxial and biaxial nematic. The symmetry consider- 
ations of those structures indicate that the leading order parameters are associated with the 
symmetry adapted base functions with angular momentum index / equal 2. That leads to the 
truncation of the expansion of the bifurcating state p\ (ft) in (2.35) after I = 2 and implies that 
the leading contribution to the bifurcation matrix d/ 2 ** (2.40) from the pair direct correlation 

(2) 

function comes from the coefficients &m,n- When the isotropic phase is taken as the reference, 
those quantities determine the bifurcation point completely. In case of r% - N B transition, p 
additionally depends on the uniaxial order parameters Aq m , which due to the non-isotropic 
structure of the reference state come into the bifurcation equation up to I = 4. 

We would like to note that both equations (2.47) and (2.58) are non-linear relations ful- 

(2) 

filled at bifurcation point for density p = p and temperature t since the coefficients c m \n and 
order parameters Aq m in general depend on t and p. The description of the method used in 
the determination of Cm]n and solutions of Eqs. (2.47), (2.53), and (2.58) are presented in Ap- 
pendix B . 

Considering the transitions from isotropic phase and using only the bifurcation equation 
(2.33), we cannot say what state is described by pi(f2) because the equations for the bifurca- 
tion point of biaxial and uniaxial states are the same (2.47). In order to obtain the symmetry 
of the bifurcating phase, we need to consult the second order equation (2.34). This approach 
will also lead to the Landau points, as is described in the following section. 

2.2.4 Exact equations for Landau point 

The Landau point, as we have indicated earlier, is the point on the phase diagram where 
four phases: isotropic, uniaxial rod-like, uniaxial disc-like, and biaxial nematic meet. Thus, 
in order to locate it, we calculate the bifurcating state, find its symmetry using Eq. (2.34), and 
determine the places where the two types of uniaxial nematic ordering, and Nu-, become 
undistinguishable [28] . The other method is based on the analysis of the duality transformation 
of the part of c 2 expansion Eq. (2.44) with angular index equal 2. The thermodynamic states 
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left invariant under this operation (self-dual) coincide with Landau points [71]. Both methods 
lead to the same conditions. We start with description of the former approach. 

As was shown by Mulder [28], in order to acquire some information about the structure of 
solution that branches of the reference state, we should calculate the eigenvectors of operator 
K 2 in (2.33)-(2.34). Let's consider transitions from isotropic phase. The bifurcation equation 
states that 

Pl (fi) = ^ Po K 2 (fl, [pj) . (2.59) 

We mentioned before, that the above has the form of an eigenproblem with K 2 being a linear 
operator acting in the space of real functions in which the orthogonal symmetry adapted set 
(C.12) was chosen as a base. In that space we can define an inner product 

(/i,/ 2 ) = J dflhifyMti), (2.60) 

with respect to which the operator K 2 Eq. (2.30) is Hermitian: 

( K 2 (Cl, [A]), f 2 ) = ( f x , K 2 {tl, [/ 2 ]) ) . (2.61) 

We also know that the eigenvalues p calculated from Eq. (2.47) are the same for the sets of 

(2) 

biaxial and uniaxial coefficients a m \n (Eq. (2.35)), which comes as a direct consequence of 
(2.37) and the trivial (isotropic) structure of the reference state. 

Now following [28] we turn to look for states that branch off the isotropic (reference) 
phase. Eigenvectors associated with the eigenvalue (2.47) are 

Xo = e Ag>(fi) + e 2 Ag(fi), 

m - o\ - (2.62) 

X2 = eoA$(n) + e 2 Ag(n). 

The e„ coefficients are chosen such that \n are normalized with respect to the inner product 
(2.60): 

( Xm, Xn ) = ^T-5m,n • (2.63) 
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They read 



where 



e 



e 2 



.(2) 
-0,2 



eg) + r 2 



eg + r 2 



(2.64) 



1 

r = - 
2 



r ( 2 ) _ r ( 2 ) 
c 0,2 c 2,2 



(2) _ 12) 
0,0 c 2,2 



4 c, 



(2) 
2,2 



(2.65) 



In the symmetry adapted base (2.50)-(2.5 1), a general solution of Eq. (2.59) is given by: 



Pl(fJ) = CtQXO + «2X2 • 



(2.66) 



We can add a normalization condition for pi by setting a 2 , + a 2 = 1 and calculate the coef- 
ficients a n using (p 2 , P re f ) = and (p 2 ,pi ) = 0. With second order bifurcation equation 
(2.34) and (2.61) we get 

Pl ( X n, K 2 (Cl, N) ) + ^pl( Xn , K 2 (fi, [Pi]f ) = . (2.67) 

Combining Eqs. (2.66) and (2.63) with the bifurcation equation (2.59) we can further simplify 
the above to obtain 

\p\ao + \po [ajj( Xo, xl ) + <£{ Xo, xl ) + 2a a 2 ( X2, Xo)] = , 

11. , (2 - 68) 

-Pi«2 + -po [«o( X2, xl ) + a§( X2, xl ) + 2a a 2 (xo, xl)\ = ■ 

The products of the eigenvectors % n can be calculated; in those equations each of them con- 
tains three symmetry-adapted functions. Using (C.19) we arrive at the following conditions: 

Pia + ^8tt 2 po f («q - a|) = , 

I (2.69) 

Pl«2 - ^87T 2 p £ «0 "2 = , 

where £ = e (e^ — 3 e 2 ). 

Now we see that for £ 7^ the allowed values of (a ,a 2 ) are: (00,02) = (±1,0), 
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(a ,a 2 ) = (±1/2, ±y/3/2), and (a ,a 2 ) 
tions forpi, Eq. (2.66), are 



±1/2, =)= \/3/2). The corresponding solu- 



(3) 
Pi 

(2) 
Pi 



±Xo , 

1 , V3 

1 v/3 
± ± ^-X2 



(2.70) 



All of the above are of uniaxial symmetry. It can be easily seen once we use the represen- 
tation of Am,n(n) in directional cosines between the laboratory and body associated tripods 
|li,i 2 ,i3| and jbi,b 2 ,b 3 j as in (2.50)-(2.51). pf' 1 does not change under rotations about 

axis 1 3 , while the p^ and p^ about 1 2 and L, respectively. The appearance of the three solu- 
tions is a natural consequence of the freedom of renumbering of the axes, i.e., the permutation 
symmetry. 

The choice of positive solutions in (2.70) for £ > entails the ordering of b 3 axes along 
1 3 (see (2.50)-(2.51)), thus marking the behaviour of molecules in rod-like uniaxial nematic 
phase, while the case of £ < 0, where negative solutions are used, stands for disk-like be- 
haviour. The one distinguished case is £ = 0, which marks the point where the difference 
between oblate and prolate nematic phases ceases to exist. In that place a second order (since 
£ = Ai = 0, as can be seen from (2.69)) transition to lower symmetric, biaxial structure 
is possible. The necessary conditions can be cast in the following form: 



J2) 
-0,2 



, for c 



(2) 
0,0 



4 2 2 > o 



c (2) 

c 2,2 



J2) 
"0,0 



2 

71 



.(2) 
"0,2 



(2.71) 
(2.72) 



The above is a non-linear equation for Landau point, relating temperature and model specific 
parameters. It proves that the location of Landau point can be found once we know the co- 
efficients Cm]n [28, 80]. If c 2 2 ^ and c 2 q / we can rewrite the above as: c^l/%o = 
1 — 2/v^ c[, 2 2 /cq 2 q. This defines a line in (cq^/c^o, c^/c^o) space that marks a boundary 
between the disc-like (when $l/%o > 1 — 2/V3 /%q) and rod-like (when c^/c^q < 

1 — 2/V3 c 2 2 /cq 2 q) states. The above analysis can be generalized further, and the existence 
of Landau points can be proven [80]. 

There exists another way of determining Landau points following from general consider- 
ations on the contribution of the terms with angular momentum index I = 2 in c 2 (f2 fi) in 



46 



Local Density Functional Theory 



the expansion (2.44). Let's write it down explicitly, up to the term with / = 2 



+c 2|0 AgJ(n' i n) + C2)2 A^(n' ~n) 



(2), 



1-1 



(2.73) 



where for clarity we have dropped the upper index (c m n = Cm, n ), and we keep in mind that 
it is an approximation, since we have truncated the expansion; however, we already know 
that it gives the leading contribution to the bifurcation. As we have mentioned, the particle 
interchange symmetry requires c 2j o = co, 2 , which can be explicitly seen once we have used 
the Cartesian representation for A; n ' n (fl) functions (see (2.50)-(2.51) and also Appendix C): 



c 2 



c o,o 



c 2 



c 



biL 



C 2 



b 3 I3 



V3i 



Co 



b 2 l 2 



sgn (c 0)0 J 



C2 



(2.74) 



where c = c 0;2 / |c , | and c 2 = c 2>2 / |c 0j0 |. 

The thermodynamic properties of the system in our approach are determined by the self- 
consistent equation (2.52), we recall that the pair direct correlation function contributes to this 
equation via term pc 2 (S7 ft), namely 



P(Q) = Z- 1 



exp 



1-1 



p / do, c r (n n)p(n ) 



(2.75) 



where p = p |c 0i o| and a reduced c 2 is denoted c r (fi 1 f2) = c 2 (f2 /| c ,o I • In order to 
obtain (2.74) and in effect the above equation, we needed to number the axes associated with 
laboratory and molecular frames. We can always renumber(permute) the axes and require that 
the phase diagram does not change. It means that p ' c r '(f2 = pc r (f2 where the 
former stands for the pc r (f2 fi) with renumbered axes. In particular, if we change the "2" 
and "3" axes of the molecules, and replace c 0) o, c and c 2 with c' , c and c 2 , where 



sgn(c 0t0 )\c 0>0 \ 
sgn(c' 0fi )cQ 

sgn(c 0fi % 



c o,o| 



6 c 2 



sgn(co,o) + 2^ + 3 c 2 , 
y/3 c [3 S5fn(c , ) - c 2 ] + 3 [s#n(c ,o) 



c 2 c 2 



[sgn(c 0fi ) + 2^3 c + 3c 2 ] (3c + v / 3c 2 ) 
-6 c + [3 s#n(c ,o) + c 2 ] 



(2.76) 



6 c + [s#n( c o,o) + 3 c 2 ] 
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then the self-consi stent equation will give the same phase behaviour of the system, whether 
we use p' , Cq , c ', c 2 or p Q , c 0j o, cq, c 2 , because p ' c r '(S7 f2) = pc r (Q f2). We have 
obtained a non-trivial duality transformation which connects the states of different density 
and temperature. The interesting points are the self-dual ones, i.e., those left unchanged by the 
above duality transformation, they fulfil c = c and c 2 = c 2 . Those are the point (c , c 2 ) = 
(-s^n(co,o)/V3, -s^n(co j0 )) and the line 

c 2 = sgn (co |0 ) - 2c / v 7 ^ . (2.77) 

This condition is the same line as the one obtained previously from bifurcation analysis and 
expressed by Eq. (2.72). It marks the cross-over between the states of prolate and oblate 
symmetry [10, 71, 81]. In the similar way we can work the other duality transformations [71]. 

The Landau points can be found by both the analysis of the bifurcating state, and the study 
of the regions invariant under the duality transformation. In the first method the Landau points 
are expressed by means of eigenvectors of the operator related to the pair direct correlation 
function, without employing any additional expansions. On the other hand the analysis of the 
self-dual points is purely geometrical and requires nothing more but the L = 2 model of c 2 and 
reasonable requirement of invariance under the permutation of the axes of the molecule-fixed 
orthogonal tripod of vectors. 

The above analysis allowed to determine the symmetry of state bifurcating from isotropic 
phase and the equation for Landau point. It is also possible to find the condition for the point 
where the character of bifurcation changes, it is presented in the next section. 



2.2.5 Equations for tricritical point 

As can be seen from Fig. 2.1, we can locate the points where the character of the bifur- 
cation changes and the corresponding first order phase transition becomes continuous. This 
behaviour is observed at tricritical point. The required condition is p 1 = p 2 = [71] which 
marks the path between schemes Fig. 2.1(b) and Fig. 2.1(c). Detailed calculations assuming 
the expansions (2.73) and (2.48) give 



3r], 



2 2 



Vl + 3— 



P 



£o - £o 2 - — sgn(c 0i0 ) e 2 - 2 (rj £ - Vo£o + -7C0) ef 



1 

P* 



~ (m 2 -ril + -c 2 )f 2 
p* 



6r)id[2r) 2 t 2 + d£\+d 



-4 772^2 



(2.78) 



+d 12 t? 2 £ 2 2 - 24 77 2 2 £ 2 2 (3 ^ £ 2 - 7/2$ ) d + 3 $ 2 - £ 2 4 )d 
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where 



(2.79) 



d-- 

p*- 1 sgn(c 00 ) - £f 
e =g [g 3 (lnZ 2:3 p*~ 1 -Cq) - g 2 (lnZ 3;3 p*~ l -c 2 
f =9 [-93{lnZ 2 ,2P*~ 1 - sgn(co fi ) + g 2 (lnZ 2 - 3 p*^ 1 - c )] , 
g' 1 =p*Z [(lnZ 2)3 p*~ l -c ) 2 - (InZ^p*' 1 - sgn(c 0fi )(lnZ 3;3 p*' 1 -c 2 
g m = lnZ m Zo t o — Zofl^m + d (2 ZnZ m , Zo,! — 2 Zo 5 i ifn + <i ( ZnZ m Zi ; i — ^i,i, m )) , 

and where all of the averages are calculated in the reference, uniaxial nematic phase, p* stands 
for bifurcation point for the case of Nu - Ng transition, Z is from Eq. (2.52) and Z mjH are 
derivatives of Z with respect to a; = (A| 2 ),X\ = (A| )> ^2 = (Aq )> x 3 = (Aq 2)' an d a l so 

= s^(c , )(A 2 mi0 ) + c (A 2 mj2 ), rj m = ^(A^q) + c 2 (A^ 2 ). 

The formulas (2.78)-(2.79) become more complicated when in the expansion of c 2 (2.44) 
terms with angular momentum index I > 2 are included and when the higher order direct 
correlation functions are taken into account, therefore, the above equations are approximate 
expressions. 



2.3 Summary 

We have presented a general formulation of the so-called Local Density Functional The- 
ory. It is a well known theoretical tool widely used for simple liquids and in the crystallization 
phenomena, applied with great success to the inhomogeneous fluids and interfaces, as well as 
in studies of the elasticity of liquid and solid crystals. It was also used in the description of 
mesophase diagrams and yielded reasonable results in comparison with simulations. 

In the present work we are only interested in the behaviour of the system close to the 
transition point. So in essence we are looking for the points where the number or structure 
of local minima of T [g] changes. One way of addressing this issue is to postulate a model 
for the direct correlation functions and use the variational principle to locate numerically the 
extrema of the intrinsic Helmholtz free energy. A first step towards this goal is to study the 
self-consistent equation (2.21) in pursuit of the points where from the stable reference state 
a new solution bifurcates. This approach allows to establish a connection between the mi- 
croscopic quantities such as molecular dimensions, interaction parameters, etc., and global, 
macroscopic properties of the system, such as the type and degree of order of the bifurcating 
state. 
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We derived exact bifurcation equations for arbitrary reference state P re f- The point of 
bifurcation was fully determined by the pair direct correlation function C2 and one-particle 
distribution function P re f, while its character (first order vs multicritical) depended on higher 
direct correlations. Then, we presented a specific form of the equations in spatially uniform 
regime for arbitrary symmetry of phase and molecules and given c 2 and P re f. The bifurcation 
equation was cast in the form of eigenproblem for bifurcation matrix uj^ of dimension not 
exceeding (21 + l)x(2l + 1), where I stands for the angular momentum index of bifurcat- 
ing state. In particular, we showed that c 2 contributes to the bifurcation point via coefficients 
of the expansion in symmetry adapted base with the same angular momentum index as the 
bifurcating representation. We also commented on the splitting of the space of solutions of 
bifurcation equations induced by c 2 , which in conjunction with with trivial structure of P re f 
leads to the same bifurcation equations for isotropic reference state. 

Then we turned to the case of molecules and phase of D 2 h symmetry. We chose the order 
parameters accordingly to the symmetry of the nematic states, and work out the bifurcation 
equations. For transitions where higher symmetric phase is taken to be isotropic, the bifur- 
cation point was fully determined by the coefficients {cm,n} of the expansion of pair direct 
correlation function with angular momentum index 1 — 2. The eigenvectors associated with 
biaxial and uniaxial phase gave the eigenequations with the same bifurcation matrix, leaving 
the symmetry of the bifurcating state to be determined by the second order bifurcation equa- 
tion. For the uniaxial - biaxial nematic bifurcation the structure of the equation was similar, 
but it contained explicitly the order parameters calculated in the uniaxial(reference) state for 
angular momentum index 2 < / < 4. In concluding sections we presented the method of find- 
ing Landau and tricritical points. The analysis of the duality transformation of the expansion 
of pair direct correlation function up to I = 2 showed the existence of self-dual regions, which 
coincided with Landau points, as was proven by calculation of the bifurcating state symmetry. 
Finally, we gave the condition for the tricritical point expressed in terms of {c m \n}- In this case 
the choice of the subspace with I = 2 in the expansion of c 2 , and the use of only this element 
from the set of direct correlation functions was an approximation, since to this equation also 
terms with higher angular index, as well as c n for n > 2 contribute. 

As was mentioned before, the real point of the first order transition lies in slightly lower 
densities (or higher temperatures) than the bifurcation point, therefore in this case, bifurcation 
diagrams are only an approximation for the true phase sequence. The true equilibrium state 
can be derived by the analysis of the minima of the free energy. This procedure, although 
numerically more complex, can give a definite answer to the question of the state that gains 
stability. 

In this work, as far as the transitions from the isotropic phase are concerned, we can expect 
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only two scenarios: to either enter uniaxial or biaxial nematic, and when the reference phase is 
taken to be uniaxial nematic only one: to stabilize the biaxial nematic phase. Those transitions 
are weakly first or second order, and while other liquid crystal mesophases are neglected, the 
bifurcation analysis is a good method for estimation of the potential possibility that phase of 
given symmetry could be realized in practice. In more complex cases a care should be taken, 
for the bifurcations do not necessary provide the correct physical results. Especially it should 
be noted once more, that certain transitions are not immediately accessible from the bifurca- 
tion analysis; only the transitions between states the symmetries of which can be connected 
by the group - subgroup relation can be straight away determined in this approach. 

Present work is based on the analysis of the self-consistent equation in the manner de- 
scribed above, that is, we focus our attention on the determination of transition densities 
and temperatures by analysing the equation for equilibrium one-particle distribution function, 
while assuming a given pair direct correlation function, which will be in fact modelled with 
the help of the mean field and virial expansion (as described in the following chapter). It is 
also possible to address the issue of determination of thermodynamic properties of the system 
by calculating direct correlation functions as, e.g., in [82]. 

In the following chapters we turn to the determination of the behaviour of the model sys- 
tems in vicinity of the phase transition mainly by finding the solutions of the bifurcation equa- 
tions (2.47) and (2.58), and also by analysing the tricritical condition (2.78) and minimising 
the free energy. To perform the actual calculations we postulate a model pair direct correlation 
function c 2 , which is equivalent to making certain approximations for the excess Helmholtz 
free energy T ex (2.22) and the reference state in the Taylor expansion (2.26). In the next chap- 
ter we start with the description of these additional assumptions, and then continue to study 
transitions to biaxial nematic for model pair potentials. 
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In the previous chapter we introduced the Local Density Functional Theory and bifurcation 
analysis with a special emphasis on the presentation of exact results. We developed the 
conditions for tricritical and Landau points following from duality transformation as well 
as from the analysis of the symmetry of bifurcating state, and derived bifurcation equations 
for transitions involving spatially uniform phases. In this chapter the solutions for those are 
presented in the form of phase diagrams at bifurcation point, for three models of the biaxial 
nematic. 

3.1 Introduction 

We start with the description of numerically tractable approximation of the excess 
Helmholtz free energy carried out through modelling the direct correlation functions, which 
is used in calculations. Then, we consider the mean field of a simple two-point interaction 
as an introduction of the basic ideas, as well as an example of the generic D 2 h - symmetric 
model giving rise to the biaxial nematic phase, Landau, and tricritical points. In this case 
we illustrate the use of the duality transformation using Eqs. (2.76) and formula (2.78). In 
following sections we consider the biaxial Gay-Berne potential and a model of the bent-core 
molecule using second order virial approximation for the excess Helmholtz free energy (2.22). 
In the first case we aim at determining the topologies of the phase diagrams as function of the 
potential parameters; we focus on the competition between molecular dimensions and various 
interaction forces in the stabilisation of the biaxial nematic phase and their influence on the 
Landau point. Then, we study models of bent-core molecules, where we seek the main factors 
that make N B stable in those systems. We use Gay-Berne interacting segments to model the 
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molecules, and include dipole-dipole interaction. Most of the diagrams presented in this chap- 
ter follow from the numerical solution of the bifurcation equations (2.47) and (2.58), where 
p is identified with the bifurcation average number density, and temperature dependence is 
stored in coefficients Cm, n and order parameters A<p m . The equations are solved in iterative 
manner, with Cm]n calculated numerically (see Appendix B). 



3.2 Models of direct correlation functions studied 

Up to now we dealt with exact equations; besides postulating the L = 2 model expansion 
(2.44) for c 2 (justified in case of bifurcations), we have made no approximations concerning 
the direct correlation functions both in the Taylor expansion (2.26) and bifurcation equations 
(2. 33), (2. 34). However, in order to find the solutions for Eqs. (2.47)-(2.58), we need a way to 
calculate c 2 (xi, x 2 ) for given x\ and x 2 . To this end we have to postulate a model function, 
and use it with orthogonality conditions for symmetry adapted functions (2.36) to calculate the 

(2) 

coefficients cJm\n in (2.44). We choose the so called low-density approximation, also known as 
second order viral, or Onsager approximation [83]. In this approach we set 

/3c 2 (xi,x 2 , [Bref]) =exp [-pv(x 1 ,x 2 )] - 1 = c(x!,x 2 ) , 

(3c n (x U ...,X n , [Qref]) =0 , U > 3 , 

where V(xi,x 2 ) stands for pair potential and c(x\, x 2 ) for second Meyer function. In terms of 
the Taylor expansion (2.26) it means setting g re f = 0, T ex [Qref] = 0, and taking T ex as 

PFex [q] = ~ J dxidx 2 g(xx) {exp [-f3V(x 1 ,x 2 )] - 1} g(x 2 ) . (3.2) 

Both choosing the state of zero density as the reference for the expansion of the excess 
Helmholtz free energy and approximating the structure of pair correlation function by low- 
density model (3.1) seem like great simplifications (see Fig. 3.1(c)). We can expect that this 
approach will give higher transition temperatures than, e.g., Monte Carlo simulations. This 
overestimation is mainly due to the lack of proper treatment of entropy. Taken all this into 
account, the phase diagrams produced using this method for spatially uniform mesophases are 
surprisingly accurate [84]. 

There exists a simple limit to the (3.2); when the molecules become very long and thin, the 
main contribution to the excess Helmholtz free energy will come from the two-body excluded 
volume. Indeed, in this case, we can assume hard interaction, i.e., V(xi,x 2 ) = oo when 
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r r r 



(a) Example of pair inter- (b) Mean field approxima- (c) Low-density approxima- 
molecular potential. tion. tion. 



Figure 3.1: Two popular approximations for the pair direct correlation function c 2 
(Figs, (b) and (c)) for an example of model two-point potential V (Fig. (a)) as function 
of intermolecular distance r (remaining degrees of freedom are held constant). Low- 
density approach in Fig. (c) models more correctly the volume which is forbidden for 
molecular centres of masses due to impenetrability of "hard core" of molecules. Also 
in comparison to mean field in Fig. (b) the contribution from the minimum of V is 
stronger. 



molecules are in contact, and V(xi, x 2 ) = otherwise, and then (3.2) is replaced by: 

P^ex[o\ = ^ J dxidx 2 q(xi)Q(£,(x 1 ,x 2 ))q(x 2 ) , (3.3) 

where is the Heaviside step function, and where £ is the contact function. Eq. (3.3) rep- 
resents well known Onsager model [83], which becomes exact for pair interactions with 
anisotropic hard core and attractive soft tail in the limit of large molecular elongations (when 
length-to-width ratio is greater than 10 : 1). 

Often even a simpler model than low-density (3.2) is used. The so-called mean field theory 
can be obtained from high temperature expansion of the Meyer function c(x 1 , x 2 ) in (3.2) 1 . In 
this theory the excess part of the Helmholtz free energy reads 

fiFex [e] = 7} J dx ± dx 2 s( x i) [/3V(x 1 ,x 2 )] q{x 2 ) . (3.4) 

For a representation of the above approximations see Fig. 3.1. 

In present chapter we will only consider the spatially uniform phases, mostly in the low- 
density approximation. As we have mentioned, the pair potential and correlation function 



'For the description of the usual method of derivation of Eq. (3.4) please refer to the Appendix E. 
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depend on the relative position and orientation of the molecules: c{x\, x 2 ) = c(fl 1 ?i — 
r 2 ). In that case the intrinsic Helmholtz free energy functional (2.26) using (3.1)-(3.2) can be 
cast as 

(3J 7 [P] = v P J dn p(n) in p(n) - W P 2 J dhdh' p(n) c (n', n)p(n) , (3.5) 



where as before p = (N)/V stands for average number density, P((t) is the one-particle 
distribution function normalized to unity, and c( Cl , f2) = J d 3 rc(n 1 £l 1 r). In equilibrium 
P(f2) fulfils the self-consistent equation (2.52), which reads 



P(fi) = Z- 1 exp 



p dn c(n,n)P(ct ) 



(3.6) 



where Z — J dtl exp p J dti c(f2 , fi)P(f2 ) ensures the normalization. In the above equa- 
tion, the temperature is present in c(fo , f2) through the formula (3.1) and p is an independent 
parameter. It is not the case when c(f2 , fi) = —f3 J d 3 r V(fl , 17, r), i.e., when mean field 
model is in force and the dependence between p and (3 is linear. Thus, the low-density approx- 
imation allows to obtain the more realistic (non-linear) phase diagrams in density-temperature 
plane, while, on the other hand, mean field works well for lattice models [85]. 

The formulas (3.5)-(3.6) can be used to calculate pressure P = pp — T jV and chemical 
potential p 



dp ■ 



pp = P - -p 2 j dn dn p(h ) c (n , n)p(n) , 

pp=i + \n P + / dnp(n)\nP(n) - p / dndn' p(n) c (n\n)p(n) . 



(3.7) 



Once we have chosen the approximation for the direct correlation function c 2 like (3.2) or 
(3.4), we can calculate the coefficients {cm'„} in (2.44), and then use Eq. (2.48) and (2.49) to 
obtain the self-consistent equations for order parameters: 



(A( 2 > ) 

\ m,nl 



dn exp 



p dn c (n,n)p(n ) 



A£!„(n) • 



(3.8) 



Then we can solve the bifurcation equations (2.47) and (2.58). 

Another possibility to obtain the behaviour of a system close to phase transition is to 
choose a trial one-particle distribution function P tria i(n) with some adjustable parameters 
and use the variational principle (2.17) for the Helmholtz free energy T [P] (3.5) to locate 
one-particle distribution function corresponding to equilibrium state by finding a set of the 
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parameters that give P tria i(Sl) which minimizes T [Ptriai]- We choose the trial function of the 
following form, inspired by Eq. (3.6), and similar to mean field (E.3): 



trial 



z- 1 



trial eX P 



5> m „Ag! B (n) 



(3.9) 



where as usual Z tr l ial 



mn ^LJn(fi) 



ensures normalization. As mentioned in 



J dSl exp a r, 

Appendix D, this is not the only model available. Employing the above formula allows to find 
the minimum of T [Ptriai] as function of a mn , but it is easier to rewrite (3.6) as the equations 
for a mn , and once those are solved, locate the minimizing set. The resulting equations have 
similar form to the Eq. (3.8), namely 



a, 



i 7T Z 



p 



dsi dsi d 3 r a£1 (Ci)c(Ci\ si, r)p 



trial 



si 



(3.10) 



is found, we can use the equilibrium one- 



eq 



Z trlal U eX P 



leg 



A 



(2) 



'SI) 



to calculate 



When the equilibrium set of parameters a r 
particle distribution function P tria i(Sl) 
the order parameters along Eq. (2.49). 

In the following section, we present the bifurcation study of a simple L = 2 model in mean 
field approximation (MFA), then, using the low-density approach we consider the models of 
biaxial Gay-Berne interaction and bent-core systems. Mostly we present the solutions of the 
equations (2.47) and (2.58) using (3.1)-(3.2) and determining coefficients Cm,„ by calculating 
numerically the integrals in (2.45). In one case we apply the model one-particle distribution 
function (3.9) and minimize the Helmholtz free energy. A more technical description of the 
way Cm)n were obtained and the methods of addressing Eq. (3.8) are presented in Appendix B. 



3.3 Mean field of L = 2 model of biaxial nematic 

Present section is devoted to the detailed analysis of the generalized Straley's [10] inter- 
molecular potential model. It is the simplest interaction giving rise to the biaxial nematic. In 
the same way as in the Maier-Saupe case, where the orientational part of the potential energy 
possesses the same symmetry as the uniaxial nematic [6], present model is D 2 h - symmetric 
and describes systems of molecules of according symmetry. It contains a complete set of four 
order parameters needed to fully describe the biaxial nematic. It allows to construct a generic 
bifurcation diagram in MFA, which is exact in this case, as function of the coupling constants, 
and also gives general ideas about the mechanisms of stabilisation of N B . Furthermore, the 
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Table 3.1: Parameters of models studied so far in relation to our model. 

Model w ,o v 0)2 v 2}2 

Straley[10] (3 (2/sqrt3)j 5 

Two-tensor [92] —1 — v^T — 3A 

Dispersion [86, 87, 88, 89, 90] -1 ±^7 -2A 2 
Amphiphilic [91] 0-1 



model exhibits tricritical behaviour, and additional analysis of its symmetries reveals the exis- 
tence of the Landau regions, as we have seen in Sec. 2.2.4 and 2.2.5. 

The generalization of the Maier-Saupe potential to the D 2 h symmetry reads [10] 

Many special cases of this model were studied previously. We will briefly describe those 
results, and indicate the relation to the present model. In the following studies, it was a priori 
assumed that the interaction is of the form (3.11) with v m>n as constants, it is the simplest 
possible approximation that can give rise to the biaxial nematic phase, more than thirty years 
after its introduction it remains a subject of studies. 

Both mean field and simulations were used in the past to deal with the Straley's interaction 
(3.1 1) or even more simplified models. Among them were dispersion [86, 87, 88, 89, 90, 91], 
two-tensor [92] and amphiphilic [91] models. Table 3.1 has the relation of those to parameters 
used here. 

For dispersion model the mean field and Monte Carlo predicted a single Landau point for 
A = 1/ "s/6, where second order isotropic-biaxial nematic transition takes place [86, 87, 88, 89, 
90]. Simulations also suggest that when v 0)2 ^ and v 2>2 = N# does not emerge. Therefore 
"the most" minimal coupling capable of producing the biaxial nematic phase is v 0,2 = and 
v 2,2 < 0. That led to the extreme simplification setting only v 2j2 = —1, for which bot MFA 
and simulations predict a direct I so - N# transition [91]. Two-tensor model was also studied 
in the context of stability of biaxial nematic and showed isolated Landau points [92]. There 
also a special case of 7 = was examined by mean field. For < A < 0.2 uniaxial-biaxial 
nematic transition was found to be of second order, while for 0.2 < A < 0.22 it was first order 
and finally for A > 0.22 a first order isotropic - N B transition took place. 

The possibility of the tricritical point on the line of Nu -N B transitions was studied for 
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two-tensor model [92] (where it was found) and on the line of I so - Ng transitions in [93]. 
Also some extension to mean field model, concerning the behaviour of order parameters and 
localizing the equilibrium state, was presented [94, 95]. 

In following sections we are taking into account a general model as described by us in [71] 
and show the bifurcation study in MFA. 



3.3.1 Bifurcation for L = 2 model 



Having chosen a model interaction of the form of (3.11), we can apply MFA described in 
Appendix E, and expressed by Eq. (3.4). In the case of L = 2 model we can present exact 
solutions for the bifurcation equations (2.47) and (2.58). The analysis which was conduced 
for L = 2 model for direct correlation function c 2 (2.44) holds for the analogical expansion 
for pair potential (3.11). We can immediately use the bifurcation equations (2.47)-(2.58) and 
conditions for Landau (Eq. (2.72)) and tricritical (Eq. (2.78)) points, keeping in mind that the 
relation between coefficients present in (2.44) and (3.11) is: c m \ n = —pv m>n . 




Figure 3.2: Mean field bifurcation diagram for L = 2 model, for sgn(vo t0 ) > 0. Tri- 
critical points are represented by thick black line and self-dual line of Landau points 
is marked with grey. Upper surface stands for bifurcation from isotropic phase, the 
lower one to Nu - Nb case. 
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Without loss of generality Eq. (3.11) can be rewritten as 



j-i 



v{ti *n) = -|u ,oN s^n(uo,o)A55(n' *n) + u A$(n' *fi) + A$(n' ^ 



/-I 



.(2), 



r-1 



(2), 



/-l 



(3.12) 



where v = v 0j2 / |i>o,o| an d u 2 = ^2,2/ I ^0,0 1 > an d due to particle interchange symmetry we 
have V2,o = i>o,2- We will use the reduced temperature t = /c_bT/(|w ,o| p)- The temperature 
for bifurcation from isotropic phase can be calculated as 



1 

t= io 



agn(v 0fi ) + v 2 + yj (sgn(v 0fi ) - v 2 ) 



Av 2 



(3.13) 



The bifurcation equation for the uniaxial - biaxial nematic transition can be cast in the follow- 
ing form: 



v 2 



be) - [sgn (u ,o) + 2 a v ] t + t 2 
sgn (t>o,o) (a 2 — be) + bt 



(3.14) 



where 



70 a =20 A 



(2) 
0,2 



15A 



(4) 
0,2 5 



70 b =14 + 20 A^ + A^ + V35A 



i 0,4 ' 



(3.15) 



70c=14-20A^ + 6A 00 



.( 4 ) 



The bifurcation diagram following from Eqs. (3.13)-(3.14) with tricritical points obtained 
using Eq. (2.78) is presented in Fig. 3.2 (in the analysis we set sgn(v 0:0 ) = 1). Fig. 3.3 shows 
the line of tricritical points (marked with thick black line on the upper surface in Fig. 3.2) in 
{vo, v 2) plane in region of rod-like states. Using the duality transformation (2.76) it can be 
projected onto the rest of the parameter space. 



3.3.2 Summary 

The L = 2 model (3.11) developed by Straley [10], consistent with the choice of order 
parameters as averages of functions Am, n (f2) for m,n = 0,2 (2.50)-(2.51), despite its sim- 
plicity provides a rich behaviour. We have shown how the bifurcation temperature changes as 
function of the potential parameters, and also solved a condition for the tricritical point. Since 
the coefficients v 0>0 , v , and v 2 can be related to the model specific, microscopic parameters, 
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Figure 3.3: v 2 versus v at tricritical point in the region of rod-like states. The self- 
dual (dashed)line v 2 = 1 — S=v o is shown as well. Tricritical temperature varies from 
t = 0.2147 for v = till t = 0.2367 at v = 0.3194, where the tricritical line meets 
the Landau line. Tick marks correspond to intermediate, equally spaced temperatures. 



the results presented here contain a doze of generality; they are valid for any expansion of the 
L = 2 type, including Eq. (2.44). Many potentials can be accurately modelled in this way, 
some of the others can be cast on the L = 2 subspace and localized on the diagram in Fig. 3.2. 
In this way one can determine the minimal, required conditions for N B to emerge. In that 
sense the Straley's model, after over thirty years, still attains its significance in the search for 
the biaxial nematic. Yet it should be kept in mind that terms of higher order in the expansion 
(3.1 1) may play an important role and the bifurcation diagram may differ from the true phase 
diagram. 

In the following part of the thesis, we consider models based on the Gay-Berne interac- 
tion: the generalized biaxial version of this potential and the system of bent-core molecules 
constructed from soft, uniaxial ellipsoids. Both cases are studied in the L = 2 model of pair 
direct correlation function (2.44), which for these potentials is a good choice. We start with 
the description of the uniaxial Gay-Berne interaction since it introduces the methodology of 
soft potentials (which is also employed in D 2 h - symmetric case) and is used later to build 
bent-core molecules. 



3.4 Biaxial Gay-Berne model 



In this section we investigate the biaxial Gay-Berne model. We begin with a definition of 
the generalized potential as was done by Berardi, Fava, and Zannoni [36], then, we present 
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the results of bifurcation analysis, and where possible compare with results published up to 
now. Since we are dealing with the Gay-Berne interaction, we start with a brief description 
of the classic version of this model [96, 97]. Next, we show the generalization to the case of 
biaxial ellipsoids, as was done in [36], and then turn to the results of our analysis. We present 
the bifurcation diagram with the results of simulations marked, and take a brief look at the 
method of minimising the Helmholtz free energy using trial one-particle distribution function 
(3.9). In the following part of the section, we examine more closely interaction parameter 
space beyond the points studied by Monte Carlo. Subsequently, in pursuit of the dominant 
factor that stabilizes the biaxial nematic phase, we analyse the influence of shape and energy 
biaxiality, then the combination of those on the bifurcation diagram, with a special focus on 
Landau (self-dual) points. 

The model has rich parameter space, and there are many possibilities of investigation of 
their contribution to the process of stabilisation of biaxial nematic phase. We have chosen one 
possible parametrization in the form of a path in the space of constants entering the poten- 
tial, which will be described in following sections. Since we are dealing only with isotropic, 
uniaxial nematic of oblate and prolate type, and biaxial nematic, the best candidates for most 
important factors influencing N B stability are located by investigating Landau points positions. 
It is so because in those places the biaxial nematic is most stable, in the sense that the region 
of stability of and Nj/_ is severely reduced in the vicinity of self-dual region while at the 
Landau point isotropic phase looses stability and system stabilizes biaxial nematic. Therefore 
significant part of the analysis is devoted to the localization of Landau points as function of 
interaction parameters, temperature, and density. Taken into account a large number of pos- 
sibilities of investigating of the interaction parameter space, we consider a couple of special, 
representative cases. 

We would like to note that we are mainly interested in the transitions to the biaxial nematic, 
the behaviour of other transitions is not studied in detail. 



3.4.1 Uniaxial Gay-Berne potential 

The main success of the so-called Gay-Berne interaction lies in the proper, effective in- 
clusion of crucial for mesophase formulation properties of the intermolecular forces. This, 
in connection with mathematical tractability and low numerical cost, made it one of the most 
intensively studied models in soft matter physics [84, 98, 99, 100, 101]. 

The uniaxial Gay-Berne interaction describes the potential energy between two cylindri- 
cally (Dooh) symmetric ellipsoids of revolution. It belongs to the class of so-called soft po- 
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tentials, which means that it includes attractive forces as well as short range, steric repul- 
sions. The key quantity, in this type of models, is a "contact distance" parameter denoted by 
cr(f2 , O, f), which depends on the orientation of both molecules, more specifically on rela- 
tive orientation 17 17, and orientation of the vector linking the molecular centres of masses 
f. For length r of the intermolecular vector r equal to cr(f2 , 17, f), the Gay-Berne energy 
is equal to 0; the orientation-dependent zeroth equipotential surface is given by a(Cl , 17, f). 
This surface, for fixed orientations of molecules chosen so 17 = 17, gives three quantities a x , 
a y , a z for the intermolecular vector direction f = (1, 0, 0), f = (0, 1, 0), and f = (0, 0, 1), 
respectively, which are identified with the length of axes of interacting ellipsoids. For these 
reasons cr(17 , 17, f) is called contact distance; it is considered to be a distance between the 
centres of masses of molecules of given orientation and relative position when the ellipsoids 
are in contact, in the sense described above. In case of D^h symmetry, when, e.g. a x = a y , 
orientation of the molecule is fully described by providing a unit vector along the axis of sym- 
metry, so we can use U! = 17 and u 2 = 17, where \ii and u 2 are unit vectors parallel to the 
axes of molecular symmetry. The contact distance reads 



a (Ui, u 2 , r) 



"5* 



[r ■ ui + r ■ u 2 J 



+ 



r ui - r u 2 l 



1 + X (Ul • U 2 ) 1 - X ■ U2 



(3.16) 



where x 



l)/(ft 2 + 1), and where k = a z /a x describes the "shape anisotropy", i.e., 



elongation of the molecules. 



The Gay-Berne potential Vqb is defined as 
Vob(ui, u 2 , r ) = e(ui, u 2 , f 




(3.17) 



where a is a constant with the dimension of length, and coefficient e(iii, u 2 , f) plays a 
similar role to the contact distance only for the strength of the interaction. It determines the 
orientation dependent depth of the minimum (as can be seen in Fig. 3.4(a)), more precisely 



e (ui, u 2 , r) = e (ui, u 2 J e •* fui, u 2 , rj 



e (Ui, u 2 , r) 



2 



(?-Ui + f ■ U 2 ) 
1 + X'( U 1 ■ U 2) 



[r ux - r u 2 j 

1 + X'( U 1 • U 2) 



(3.18) 
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■e - 




"face-to-face" 



2r n 



3r n 



4r n 




-2-1012 



(a) Gay-Berne potential Vgb as a function of length 
r of the vector linking centres of masses r = rr , for 
r > cr(ui,U2,f) — <7o and for fixed f. The depth 
of the well e = e(ui, U2, f) and the position of the 
zero of Vgb at r = ro = tr(ui, u 2 , f) in face-to-face 
configuration are marked. 



(b) Equipotential surfaces for molecu- 
lar elongations 5 : 1, v = 1, ft = 2, 
£x/cz = 4. Lines indicate Vgb = 0.0 
(most inner one), and counting from 
most outer one Vgb = —0.1, —0.5, 
-0.9, -1.3, -1.7, -2.1, -2.5. 



Figure 3.4: Uniaxial Gay-Berne interaction Fig. (a) and surfaces of constant potential 
for parallel molecules (ui = u 2 ) Fig. (b). 



where 

e(ui, u 2 ) = e [1 - x 2 (ui -u 2 ) 2 ] 2 , 




(3.19) 



and where jl and v are empirically chosen exponents. The constants e x = e y and e z are pro- 
portional to the depth of the minimum of Vgb(ui, u 2 , f ) when u x = u 2 and when f = (1,0,0) 
and f = (0, 0, 1) respectively. In that sense they describe the strength of the interactions in the 
directions of main axes of the molecules. 

Fig. 3.4(a) presents an example of the Gay-Berne potential. It contains the attractive "tail" 
which behaves like 1/r 6 and models the attractive interactions between the molecules, like Van 
der Waals or dispersion forces. Those have longer range than the repulsive part of the potential 
~ 1/r 12 which models the impenetrable "core", i.e., the excluded volume that is forbidden 
for the centres of masses of molecules due to the repulsions on the shorter distances of, e.g., 
Coulomb origin. The most important feature of the Gay-Berne potential is the introduction in 
the consistent way of the contact distance and analogical parameter of the interaction strength. 
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They both depend on the relative orientation and position of the molecules; the former induces 
a shape dependent zeroth equipotential surface, while the latter introduces anisotropy on the 
attractive forces. In other words, the Gay-Berne potential opens a possibility of studies on soft 
ellipsoids of arbitrary dimensions and interaction strengths given in directions of main axes. 

The quantities in Fig. 3.4(a) have some dimension. In the analysis we need to operate 
with dimensionless objects, e.g., in Eq. (3.2) the function under the exponent does not have a 
dimension. It is also important from the numerical point of view. A natural set of units in case 
of Gay-Berne potential is associated with constants o"o and e present in (3.16) and (3.19). We 
employ a and e as the unit of distance and energy, respectively. The set of dimensionless 
quantities that will be used reads 

Vqb = V GB /e , r* = r/a , p* = pa%, 
t* = k B T/ e , /f = /V y/e a , 

where fi is the magnitude of a dipole in a dipole-dipole interaction which will be used later. 
In case of uniaxial Gay-Berne a = y/2a x . For the description of the way the dimensions are 
treated in the numerical integration procedure please consult Appendix B . 

We should note that Vqb (3.17) has another zero, apart from the one at r = cr(ui, u 2 , f) 
visible in Fig. 3.4(a). It is located atr = r — 2a and is considered to be "unphysical" together 
with whole branch of Vqb f° r r < 0"(ui, U2> r) — a . This part of the potential is separated from 
the one presented in Fig. 3.4(a) by an infinite barrier, and has no known physical application. 

On the historical account, the formula (3.16) was known as early as in 1972, originally 
developed by Berne and Pechukas [96]. They considered a "Gaussian overlap model" which 
could be understood as mathematically tractable method allowing for short range, repulsive, 
dependent on shape of molecule forces to be calculated. This is achieved by localizing on each 
molecule a three dimensional Gauss function, which can be seen as a mass distribution, and 
appropriately integrating to get er(ui, u 2 , f). Then, in 1981, the Lennard- Jones model [102] 
was modified by Gay and Berne to include er(ui, u 2 , ?) and e(u 1; u 2 , f) as in (3.17) [97], and 
ever since was known as the Gay-Berne potential. In mid-90s it was extended to the molecules 
of different elongations [103] and generalized to the biaxial, D 2 h - symmetric ellipsoids [36]. 
Now we turn to the description of the latter, the Gay-Berne potential between two molecules 
with three different axes. 
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3.4.2 Biaxial Gay-Berne interaction 



The biaxial version of the potential (3.17) has the same form, namely [35, 36] 



r - cr(f2i, f2 2 ,f) + cr c 

o_c 

r - cr(Qi, Cl 2 ,r) + a c 



12 



(3.21) 



where e and cr c are constants. Clearly, now the interaction depends on the complete orien- 
tation of the frames associated with the molecules; we need to supply three Euler angles in 
order to define the state of the molecule. Parameters a(fli, fl 2 , ?) and e(f2i, fl 2 , r) have ex- 
actly the same interpretation as in the uniaxial case. In order to define the contact distance in a 
general manner, we need to introduce a so-called overlap matrix A, the eigenvalues of which 
are proportional to er(f2i, £l 2 , f). It is built with the help of shape matrix S, which contains 
the molecular parameters, namely S ab = 5 ab a a , and reads 



A(S7i, n 2 ) = M T (fii)S 2 M(fii) + M T (ft 2 )S 2 M(S7 2 ) 



(3.22) 



where M(f2) stands for the matrix of the operator of rotation from laboratory to molecular 
frame. So, the overlap matrix A is built by rotating the squared shape matrix to the frames 
associated with the molecules, in this way the dependence on the molecular relative orienta- 
tions is introduced. The dependence on the relative position f is incorporated by calculating 
the eigenvalue of A -1 for the vector f , which leads to the contact distance 



'1/2 



(3.23) 



Similarly as in the case of the uniaxial Gay-Berne, the zeroth equipotential surface is defined 
by cr(f2i, f2 2 , f ), from the above equation. As before, the axes of the molecules a x , a y , and a z 
axe given by 



(3.24) 



where % = x,y, z, r x = (1, 0, 0), r y = (0, 1, 0), r z = (0, 0, 1). When both molecules possess 
the same orientation, the three relative configurations of the molecules with f equal to f x = 
(1, 0, 0), i y = (0, 1, 0), and i z = (0, 0, 1) will be called side-to-side, face-to-face, and end- 
to-end, respectively. Depending on the energy coefficient e(f2i, Cl 2 , f) one of them gives the 
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global minimum of the potential U (3.21) (see Fig. 3.5). The equipotential surfaces in planes 
corresponding to these configurations are presented in Fig. 3.6. 

The interaction strength coefficient is defined in the analogical way as in the uniaxial case: 



e(f2i, f2 2 , f) 

e(ni,n 2 ) 

e(fii, f2 2 , f) 



(3.25) 



2r r B~ 1 (ni,n 2 )r 



where the matrix B contains parameters of the potential strength (energy is also "biaxial") and 
is defined in a manner similar to the overlap matrix [36]: 



B(fii, n 3 ) = M r (J7 1 )EM(^ 1 ) + M T (ft 2 )EM(fJ 2 



(3.26) 



where E ab = § ab (eo/ea) 1 ^, and v, jx are empirically chosen parameters. 

The energy and length units are set by constants e and a , respectively (see (3.20)). Also 
following Zannoni [35], we set v = 3, fi = 1 [104], and a c = a y . It should be noted that 
the potential (3.21) for D^h - symmetric molecule when, e.g., o x = a y , reduces to uniaxial 
Gay-Berne (3.17) when a c = a = \plo x . 
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Figure 3.5: Dimensionless biaxial Gay-Berne potential U* = U/e against re- 
duced r* = r/a for (a x , a y , a z ) = (1.4, 0.714, 3.0), and two sets of interaction 
strength parameters. On the left the side-to-side configuration minimum is deep- 



est, \(^xi £yi ^-z, 



7-7, 1.0, 0.2), on the right face-to-face attraction is strongest; 



7.0, 1.4, 0.2). 
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-2-1 1 2-3-2-1 1 2 3-2-1 1 2 

Figure 3.6: Equipotential surfaces for the biaxial Gay-Berne potential, in three 
planes for interaction parameters (a x ,a y ,a z ) = (1.4,0.714,3.0), (e x ,€ y ,e z ) = 
(1.0,1.4,0.2). 



The interaction U in (3.21) has three shape-related, and three energy-related parameters, 
those are respectively (o x ,a y ,o z ) and (e x ,e y ,e z ). They contribute to the two coefficients 
giving rise to two sources of biaxiality: e(fii, Q2, ?) and cr(f2i, f2 2 , ?). The latter is related to 
the position of the zero of the potential and is associated with the dimensions of the molecules. 
The former gives the anisotropic depth of the minimum of the potential and corresponds to the 
biaxiality of the attractive forces. Less precisely we can say that the shape related parameters 
act along the horizontal axis in Fig. 3.5, while energy related parameters act along vertical one. 
Despite the fact that (e x , e y , e z ) and (cr x , a y , a z ) are both just constants entering the potential, 
we will call them energy parameters and shape parameters, respectively. We should keep in 
mind that they enter the excess Helmholtz free energy (3.2) via pair potential V(xi, X2) in non- 
linear manner, thus making the straight-forward interperetation of (e x , e y , e z ) and (a x , a y , a z ) 
non-trivial. 

It is useful to introduce a biaxiality parameter which can be treated as a measure of how 
much a given set of {0^} or {e.j} differs from the uniaxial, - symmetric case. There are 
several such quantities known in the literature. We have chosen to use the one following from 
the surface tensor model [105, 106] developed in [107]. It is based on the assumption that the 
mean field effective potential V e ff(x) (see Appendix E, Eq. (E.3)), which gives the equilibrium 
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one-particle distribution function P eq (x) = Z^ 1 exp [—/3V e ff(x)], can be written as a sum of 
mean torque potential V ext and internal potential. The latter accounts for the variations of the 
molecular interactions due to bond rotations, while V ext describes the external interactions 
acting on the given molecule in the sample, giving rise to the nematic ordering. Now the 
mean torque potential is expanded in the base of Wigner matrices Dm,n(&) [73, 79] (for brief 
description of these functions see Appendix C.l), following [107] 

2 

V ext (fl) ~ T 2 > n * D^h) , (3.27) 

n=-2 

where the T 2 n are the components of the interaction tensor, which is called a surface tensor. 
They read 

T ^ = - [ dSDjfhn), (3.28) 



where S denotes the surface of the molecule. Above equation is a consequence of the as- 
sumption that each surface element has a tendency to align parallel to the director in Nj/+ and 
perpendicular to it in N^_ phase. This means that V ext (Cl) ~ f s dS |(3 cos 2 (3 — 1), which is a 
base for Eq. (3.27). Coefficients T 2,n in general depend on the order parameters and molecular 
anisotropy. The molecular biaxiality parameter is defined as 

A = t 2 ' 2 /T 2 ' . (3.29) 

In case of hard parallelepiped with dimensions a, b and c, along axes x, y, and z, the above 
can be expressed as [107] 

A = ( 3 /2) 1/2 ,° (a ~ 6) r . (3.30) 
v ' J c (a + b) - lab 

This parameter will be used in this study (in place of a, b, c we will use axes of ellipsoids). It 
has some useful properties. Firstly, assuming that a / 0, 6 ^ 0, c ^ 0, the only case when 
A = occurs for - symmetric shape when a = b, for other uniaxial shapes when a = c 
and b = c, A takes values of \/3/2 and — a/3/2, respectively. When c becomes very large, 
A behaves like (a — b) / (a + b), which means that it does not become zero in this limit. It 
reaches vanishingly small value for c tending to zero. 

We choose the biaxiality parameters A CT and A e respectively for shape and energy, along 
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(3.30), they read 

K = (3/2) 



1/2 °z {<?x - CT y ) 

,1/2 e z (£g ~ 



(3.31) 

■ 

A e = (3/2) J 



i^x €y) 26^6^ 

For hard biaxial ellipsoids there is a simple condition for self-dual (Landau) point, it is the 
mentioned earlier square root rule, obtained from Eq. (2.72) [29, 80]: 

°x = y/CTyVz ■ (3.32) 

It defines a line in the plane of axes ratios, e.g., [a x j a z , a x / a y ) (see Fig. 3.7(d)), and along this 
line we can represent \ e and \ a as functions of e x /e y and a x /a y , respectively (see Fig. 3.7(c)), 
namely: 

yd _ (3/2) l/2 O x ja y - 1 _ ^ ^ 

o'x/o'y — 2 (a x /(Ty) +1 

K d - (3/2) 1/2 - ^TVi ; • ( 3 - 34 ) 

e x /e y -2(e x /e y ) +1 



In order to help to visualize the way biaxiality parameters (3.31) work, we present Figs. 
3.7(a)-(b). They show one possible method of changing X a (the same holds for \ e and {ej}) by 
varying the ratios of axes o x jo z and cr x /a y . Similar method was used in present work. We start 
in the region of biaxial rod-like shapes when a x < a y <§; a z (Fig. 3.7(a)). Then, by increasing 
the ratio a x /a y for given a x /a z (moving to the right on Fig. 3.7(a)) we decrease \ a to reach 
uniaxial, prolate molecule for o x = a y (A CT = 0, white colour in Fig. 3.7(a)). Continuing 
to the disc-like region, where a z > a x ^> a y (Fig. 3.7(b)), we start with the biaxial oblate 
molecule (dark areas), and by increasing the ratio o x /a z for given a x /a y (moving upwards on 
Fig. 3.7(b)) we reach uniaxial, disc-like molecule for a x = a z and \ a = a/3/2 (white areas). 
As we can see, in rod-like regime, vanishing A CT indicates uniaxial molecule, while disc-like, 



Dooh - symmetric shape entails \ a = y 3/2. 

The biaxiality parameters (3.31) for given shape and energy parameters {<t,} and {e^} are 
not defined in a unique way. In order to use A e and \ a to parametrize the bifurcation diagrams, 
we need to make certain assumptions concerning {a;} and {e^}, and choose the way of varying 
them. Assuming that for the biaxial Gay-Berne, there exists a self-dual geometry analogical 
to the square root rule (3.32) obtained for hard biaxial ellipsoids [29, 80], we can expect a two 
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CT X /(J y 

(a) Projection of shape biaxiality parameter X a 
to the plane of axes ratios <j x /<j z and a x /a y in 
the rod-like regime. Darker regions indicate 
lower (negative) values of A CT , (dashed)lines of 
constant \ a are drawn and labelled. The uni- 
axial limit at \ a = for <j x /<j v = 1 is marked 
with white. 




2.0 3.0 4.0 



a x /c y 




2.0 2.5 3.0 3.5 4.0 



CT X /Cy 

(b) Projection of shape biaxiality parameter X a 
to the plane of axes ratios a x ja z and <j x /a y 
in the disc-like regime. Whiter regions indicate 
larger values of X a closer to the uniaxial shape at 
Act = \/3/2, for <7 x ju y = 1. Dashed lines with 
labels indicate the regions of constant shape bi- 
axiality parameter. 
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(c) Biaxiality parameter at Landau point for (d) Condition for Landau point for hard ellip- 
hard ellipsoids (at self-dual geometry) X^ d as soids a x = ^Jo y a z [29, 80]. 
defined by (3.33). 



Figure 3.7: Closer look at the shape biaxiality parameter \ a (3.31). On the top the 
density plots of A CT in oblate and prolate regimes (on the right and left, respectively) 
as function of axes ratios. Below on the right the dependence of axes ratios at Landau 
point for hard ellipsoids, on the left \ a along this curve. 



dimensional surface of Landau points parametrized by A* and Af (3.33)-(3.34). The deter- 
mination of this surface is interesting, but it remains to be done in future studies. Presently, 
we investigate how, by varying the biaxiality of the molecule and forces, we intersect the 
self-dual plane. In other words, we choose a path in six dimensional space of {<7j} and {e^}, 
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p* p* 

Figure 3.8: Two bifurcation diagrams for parameter sets studied in [35] for \ a = 
0.58 ((cr x , a y , a z ) = (1.4, 0.714, 3.0)). On the left the model where side-to-side 
configuration is preferred A e = —0.06 ((e^, e y , e z ) = (1.7, 1.0, 0.2)), on the right 
the face-to-face minimum is deepest; \ e = 0.042 ((e x , e y , e z ) = (1.0, 1.4, 0.2)). 
Transition points discovered by Monte Carlo in [35] are marked. 



in the manner described above, and along it we calculate the density-temperature bifurcation 
diagram. Due to chosen parametrization, we cross (Af d , X s a d ) plane in a single point. Since 
the formulas (3.31)-(3.34) follow from simple models, we cannot assume a priori that they 
possess any physical meaning, especially taken into account the complicated manner in which 
{crjjand {ei} enter the potential in present model. However, we can investigate if (Af d , X s a d ) is 
in any way related to the surface of Landau points obtained from (2.72) for biaxial Gay-Berne 
interaction. 

The sign of both parameters A CT and A £ is of importance, for it reflects the ratio of a x / a y and 
e x /e y . Let's imagine two identical, rod-like, biaxial ellipsoids with a z > a x > a y (A CT > 0). It 
seems natural to expect that in the face-to-face configuration the interaction U will be strongest 
(see 3.5), which means that e x > e y > e z and in turn implies A e > 0. Monte Carlo study of 
the potential (3.21) did not predict a stable biaxial nematic phase in that case [35], instead the 
simulations showed N# for potential favouring the side-to-side configuration, namely when 
\ t < while \ a > 0. The difference between these two situations is depicted in Fig. 3.5, 
where physical (for r > a(fli, Q2, £ ) — cr c ) regions of U are shown. In present study we refer 
to the case where e x > e y > e z (A e < 0) and a x > a y (\ a > 0) by "strong lateral interactions 
model". 

The potential of the form of (3.21) can be immediately put in place of V into Eqs. (3.1)- 
(3.2) and used in calculations as described in Appendix B. Then, bifurcation equations (2.47), 
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(2.58) can be solved, and the bifurcation density p* = pi and temperature t* (see (3.20)) can 
be determined, provided we had calculated coefficients Cm,n and the reference state. Unless in- 
dicated otherwise, in the following p* will stand for the dimensionless density, following from 
the solutions of the bifurcation equations. We start with the comparison of the bifurcation 
results with those obtained from Monte Carlo in [35]. 

The interaction (3.21) was studied by Berardi and Zannoni by means of Monte Carlo sim- 
ulations in the isothermal isobaric ensemble for 8192 soft, biaxial ellipsoids. They took into 
account one set of shape parameters in the rod-like regime: (a x , a y , a z ) = (1.4, 0.714, 3.0) 
(A CT = 0.58), and following sets of potential strength coefficients (e x , e y , e z ): (1.0, 1.2,0.2) 
(A e = 0.025), (1.0,1.4,0.2) (A e = 0.042), (1.4,1.0,0.2) (A e = -0.042), (1.7,1.0,0.2) 
(A e = —0.06). Only for the last one the stable biaxial nematic was found. We investigate 
the issue of A CT and X e of opposite signs in later subsections, presently we show the compar- 
ison of our results with those obtained from simulations. In Fig. 3.8 we plot the bifurcation 
diagrams for parameters studied in simulations [35] and mark the transition points discovered 
there. It should be noted that the Monte Carlo study [35] did not include the long-range cor- 
rections, i.e., the potential involved a cut-off radius. For this reason the comparison is only 
qualitative, and the bifurcation parameters following from low-density differ from simulations 
results, especially for I so - N;y phase trasition. Once the effects of cut-off are compensaed 
in Monte Carlo, this comparison is improved [108]. Table 3.2 shows the critical densities and 
temperatures from simulations in comparison with bifurcation for I so - N[/ + and Nj/ + - Nb 
transitions. 

3.4.3 Contribution from precise reference state 

In Sec. 3.2 we have pointed to the minimisation of the Helmholtz free energy T [P] (3.5) 
with help of the trial one-particle distribution function (3.9) as an alternative method of obtain- 
ing a phase diagram. The above figures show how Monte Carlo results differ from bifurcation, 
which, besides the mentioned long-range corrections, may be caused by the L = 2 model 
of pair direct correlation function c 2 (2.44) used in the process of obtaining of the reference 
one-particle distribution function P ref along Eqs. (2.53)-(2.55), or by the loss of accuracy in 
the bifurcation analysis. The minimisation technique addresses these issues. In this method 
we do not expand c 2 and study the minima of T [P] directly. Therefore it is worthwhile to 
present the comparison of those two approaches. We do so in present section. 

The phase diagram presented here follows from postulating the trial one-particle distribu- 
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Table 3.2: Transition points obtained in Monte Carlo simulations (MC) by Berardi 
and Zannoni (model parameters as on the left plot in Fig. 3.8) compared to the results 
of bifurcation analysis in low-density approximation. 





density p* 


temperature (MC) 


bifurcation temperature t* 


Iso - 


N u+ 0.271 


3.20 


6.97 


Nc/+ - 


- Nb 0.291 


2.90 


3.50 



tion function Eq. (3.9) and finding the set of parameters a m n that minimises the Helmholtz 
free energy (3.5), by solving the self-consistent equations (3.10), and locating in this way the 
equilibrium state. Once all such points are found, we can easily decide which one of them 
is a global minimum. Having determined the set of parameters a m ^ n giving equilibrium one- 
particle distribution function, we can calculate the order parameters using the definition (2.49) 
and derive the phase transition points. It is also possible to find the coexistence (or more pre- 
cisely a mechanical stability) regions by looking for points where the pressures and chemical 
potentials are equal. We have done that analysis and the results for Njj - Nb transition can be 
seen in Fig. 3.9, where the phase transition as well as bifurcation lines are shown. 

Since, by the investigation of the order parameters, the uniaxial - biaxial nematic transi- 




0.16 0.20 0.24 0.28 



p* 

Figure 3.9: Phase diagram for r% - N B transition following from minimisation of the 
Helmholtz free energy (3.5) by solving of the equations (3.10) (solid line), and from 
bifurcation analysis (dashed line) for potential parameters studied in [35]: A CT = 0.58, 
A e = -0.06 ((a x , a y , a z ) = (1.4, 0.714, 3.0), (e x , e y , e z ) = (1.7, 1.0, 0.2)). 
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tion has been found to be second order and because in this case the bifurcation point is a point 
of phase transition, the figure presents the difference between the bifurcation point obtained 
with the help of reference state P re f calculated using of L = 2 model of C2 (2.44) and the 
exactly obtained P re f, as described above. In other words, the difference between solid and 
dashed line in the Fig. 3.9 is the estimate of the influence of all coefficients {cm,n} with I > 2 
in P re f (2.55) on the phase transition point obtained from bifurcation equation (2.33). The 
comparison of bifurcation and phase diagram indicates that the critical values of density and 
temperature obtained by those two methods are complementary. This proves that the leading 
contribution to the transition density and temperature comes from terms in the expansion of c 2 
with angular momentum index equal 2, and that, at least, for transitions studied here, the bi- 
furcation analysis gives reasonable values of critical parameters. The numerically much more 
expensive free energy minimisation shifts the diagram, but does not change the overall phase 
behaviour. 
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Figure 3.10: Path of varying of shape biaxiality parameter A CT (thick, solid line drawn 
with gradient) in plane of axes ratios (a x /a y , a x /a z ), and the projection of A CT . Lan- 
dau point resulting from present analysis (shown in Fig. 3. 17(a)) is marked with circle 
at (2.57, 0.44), the self-dual, Landau region for hard ellipsoids (a x = ^/a y a z ) is rep- 
resented by dotted line. Uniaxial regimes for X a = 0.0 and X a = a/3/2 are marked 
with white. Thin, solid lines indicate regions of constant X a . 

As can be seen from Figs. 3.8 and 3.9, the dependence of temperature t* on density p* at 
bifurcation is almost linear; the variation of density does not change the topology of the phase 
diagram. Therefore we can choose some value of p* (keeping in mind that assumptions made 
in Chapter 2 for the Taylor expansion (2.22) may prove incorrect if it is too high) and pursuit 
the effects of biaxialities on bifurcation temperature at fixed density. That will be the subject 
of the following section. 
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3.4.4 Exploring the effects of biaxiality 

In the biaxial Gay-Berne potential there are two sources of biaxiality identified with the 
anisotropic forces strengths and the shape of molecules. The deviation from uniaxial symme- 
try of those is expressed by parameters A £ and X a , respectively. Currently we will show how A e 
and \ a and combination of these change the bifurcation diagram. We present the dependence 
of dimensionless temperature t* on dimensionless density p* at bifurcation point for fixed A e 
and A CT , as well as t* as function of A CT or \ € for fixed p*. In the latter case the diagrams are 
calculated for p* = 0.18, unless indicated otherwise. It was chosen to give packing fraction, 



phases: 0.4 < rj < 0.6). 

As we mentioned before, in order to parametrize the bifurcation diagrams by A e and A CT , 
we need to make additional assumptions concerning {a.;} and {e^}. We have decided to hold 

(Tj = const and = const (which set the density and temperature scale), and q > 0, 
(Tj > 0. The sets of energy and shape parameters we refer to are listed in Table 3.3, they 
are denoted by (a) e -(i) e and {a) a -{h) a , respectively. As a starting point we have chosen the 
set studied by Berardi and Zannoni [35]: (e) £ , (b) a 2 . That fixed £\ °i an d J2i e i- From this 
point, in order to decrease the shape biaxiality in a consistent manner, we choose to decrease 
the difference between a x and a y while holding o z constant, reaching uniaxial shape (a) a for 
\ a = 0.0. In other direction we increase o z — o y while keeping a x = const up to a moment 
where Landau point is found at (g) a . From this set we make the molecule more disc-like, by 
making a x — a z and a y smaller, to reach uniaxial, oblate ellipsoid (h) a for A CT = \/3/2. This 
path, in the plane of axes ratios (a x /a z , a x /a y ), is depicted in Fig. 3.10. As we can see we 
are crossing the self-dual line for hard D 2 h - symmetric molecules once. The Landau point 
following from present analysis is also depicted. 

In case of {e^} we do not cross the boundary between the "disc-like" and "rod-like" en- 
ergy parameters sets. We only follow the path linking models of strong lateral interactions and 
cases where face-to-face attractions are strongest. That is, starting from (d) e , being the point 
of uniaxial energy where A e = 0.0, and holding e z = const, we increase energy biaxiality 
by increasing \e x — e y \ to reach (a) e and (i) e . Along this path we cross the point studied in 
simulations [35] (e) e , and also the model fulfilling square root rule (3.32) for {e^} (h) e . With 
the above sets of shape parameters, side-to-side and face-to-face configuration give deepest 
minimum of the potential when e x > e y and e x < e y , respectively. 

The above method of changing A e and A CT was used in every case. 

2 We could have chosen the other one with (e x , e y , e z ) = ( 1.7, 1.0, 0.2), the only difference between those 
two is relatively stronger interaction in the latter case. 
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Table 3.3: Biaxiality parameters A e , \ a and interaction constants sets {ej}, {(J,}. 



parameter set 


^■xi ^yi ^z 


K 


parameter set 


O ' z 


<v 


We 


1 o n ^ n o 
i.y, u.o, u.z 


— U.Z410 


{a)a 


l.UO ( , l.UO I , o.U 


n n 
U.U 




1.84,0.56,0.2 


-0.1974 




1.4, 0.714, 3.0 b2 


0.58 




1.4,1.0,0.2 


-0.042 


(c)„ 


1.4,0.637,3.077 


0.6412* 


(d) e 


1.2,1.2,0.2 


0.0 


(d)„ 


1.4,0.614,3.1 


0.6596 


(e)e 


1.0,1.4,0.2 bz 


0.042 


(e)a 


1.4,0.574,3.14 


0.6919 


(/)« 


0.9,1.5,0.2 


0.0662 


(/), 


1.4,0.564,3.15 


0.70 


(9)e 


0.8,1.6,0.2 


0.0942 


In). 


1.4,0.544,3.17 


0.7163 


(h)e 


0.6,1.8,0.2 


0.1750** 




2.35,0.414,2.35 


^3/2 


«« 


0.5,1.9,0.2 


0.2415 









(*) self-dual point for hard biaxial ellipsoids (a x = ^a y a z ). 



(**) fulfilling square root rule for energy parameters (e x = v /e^eT). 
(bz) parameters studied in Monte Carlo simulations in [35]. 



If we fix X e = for (d) e , we can expect that the system of Z^ooft - symmetric, rod-like 
molecules with axes (a) CT (A CT = 0.0) stabilizes N[/ + at the expense of isotropic state, while 
for uniaxial disc-like ellipsoids, when \ a = a/3/2 for (h) a , most probably I so - N[/_ tran- 
sition occurs first. Thus somewhere along the path of varying molecular shapes, described 
above, leading from \ a = to A CT = a/3/2, we should cross the boundary between oblate 
and prolate nematic states, and hence stabilize the biaxial nematic. What happens at the point 
where we cannot tell whether the nematic phase is of disc-like or rod-like type, and how this 
point is realized, and how bifurcation temperature and density change along the chosen path 
in {ei},{(Ti} space, will be discussed in the subsequent part of the section. 

There are some matters which are of interest, and we take them systematically into ac- 
count. Firstly, there is the question of opposite signs of A CT and A e , which was pointed out by 
Berardi and Zannoni in [35]. Next, we can wonder how the change in the shape or energy 
biaxiality separately vary the transition temperatures. Then, we can ask ourselves how those 
two parameters concur in the creation of biaxial nematic, and whether it is possible to induce a 
Landau point in this model by changing A e and X a , and if so, then how the position of that point 
changes with the biaxialities. Since we know the relation between axes of biaxial ellipsoids 
at self-dual point for the models of purely repulsive forces (3.32) (hard molecules [29, 80]), 
once we find this point in present model, we can estimate the influence of attractive (soft) 
interaction on its position, but more importantly we can check whether a similar condition is 
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(a) Temperature versus density at N[/+ - Nb 
bifurcation, for some shape biaxialities. 



(b) Bifurcation temperature as function of shape bi- 
axiality \ a . I so - Njj + and - Nb branches are 
plotted at p* = 0.18. 



Figure 3.11: Bifurcation diagrams for energy biaxiality A e = 0.0 ((d) e in Table 3.3). 
X a varies along thick, gradient path in Fig. 3.10. 



in force for parameters {cr^} and {ej}. 



Molecular shape effects 



Firstly, we take a look at the influence of shape biaxiality on the bifurcation temperature 
for the rod-like uniaxial to biaxial nematic transition. We change A CT , while keeping {e^} 
fixed, along the path depicted in Fig. 3.10 (thick, gradient line). Figure 3.11(a) represents 
bifurcation lines in (p*,t*) plane for different values of A CT , and since we are interested in 
shape dependence alone, the energy biaxiality A e = 0.0 (for (d) e in Table 3.3). As expected, 
the increase in the molecular biaxiality leads to higher temperatures for Nu+ - Nb transition, 
chances for observing the biaxial nematic phase are increasing. 

Fig. 3.1 1(b) shows the bifurcation diagram, where temperatures are plotted as function of 
shape biaxiality both for transition from isotropic phase and for Nu+ - Nb case. We can see 
that the line of uniaxial - biaxial nematic bifurcation gets closer to the I so - N u+ one, with a 
prospect of meeting it eventually. 
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(a) Bifurcation temperature versus density for (b) Bifurcation temperature as function of en- 
Nf7 + - Nb transition for a couple of energy bi- ergy biaxiality. I so - Njj + and Nj/ + - 
axialities, transitions are plotted for p* = 0.18. 



Figure 3.12: Bifurcation diagrams for uniaxial molecule shape; X a = 0.0 ((a) CT in 
Table 3.3). A e is changed, as described before, by varying \e x — e y \, while e z = 0.2. 



Biaxiality of soft part of the potential 

Now we show how the biaxiality stored in (e x , e y , e z ) affects the transition point. For fixed 
\ a , we change A e by increasing \e x — e y \ while keeping e z = 0.2. In Fig. 3.12(a) we present 
bifurcation lines (t* as function of p*) for N u+ - N B transition, plotted for different values of 
A e and for uniaxial molecule (A CT = 0.0 for (a) a in Table 3.3). As in the case of Fig. 3.1 1(a), 
the dependence of the temperature on density in Fig. 3.12(a) is almost trivial (linear), and is 
also weaker than in Fig. 3.1 1(a). 

The Figure 3.12(b) shows the bifurcation diagram with isotropic, N u+ , and N B regions 
as function of X e , for uniaxial shape A CT = 0.0 for (a) CT . It can be seen that for low shape 
biaxialities it is required to supply increasingly high energy biaxiality in order to produce N#. 
As was mentioned earlier, X e was changed by increasing the difference between e x and e y , and 
essentially we did not cross the boundary between "disc-like" and "rod-like potential" regimes, 
that is the reason why in Fig. 3.12(b) only prolate nematic phase is observed. Nevertheless, 
the plot suggests that it is possible for the lines of I so - N[/ + and Nu+ - Nb transitions to 
meet somewhere, for higher energy biaxiality. 

From Figures 3.11 and 3.12, we can see that the main factors affecting the diagram are 
the biaxiality of molecule A^ and the biaxiality of potential A e . Thus we can wonder how the 
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Figure 3.13: Bifurcation diagram for N;y + - N# transition for negative and positive 
energy biaxiality A e , (e x , e y , e z ) = (1.77, 0.63, 0.2) and (e x , e y , e z ) = (0.63, 1.77, 0.2), 
respectively. \ a increases along the path in Fig. 3.10 (thick, gradient line), p* = 0.18. 



combined shape and energy biaxialities influence the bifurcation temperature. 



Concurring biaxialities 

Turning our attention to the question of the influence of both biaxialities, now we allow A e 
and A CT to change simultaneously. It is done as before: shape biaxiality varies along the path in 
Fig. 3.10, and A e by increasing \e x — e y \ for e z = 0.2. Firstly we consider the case of positive 
and negative A £ , while keeping A CT > 0. Figure 3.13 shows the bifurcation temperatures as 
function of shape biaxiality for two cases of A e = —0.16 and \ e = 0.16, for ]%_,_ - N B transi- 
tion. This plot shows the difference between opposite signs of A CT and A e . It can be seen that 
the branch corresponding to negative energy biaxiality is above the positive one almost every- 
where, except the point where A CT = 0.0 (Fig. 3.12(b) is symmetric under the change of sign 
of A e ), and the region where shape biaxiality becomes (relatively to A e ) high enough to take 
control over the transition. A more detailed study on the concurring biaxialities is presented in 
Figs. 3.14(a)-(b), where the lines of bifurcation for N[/ + - Nb transition are plotted as function 
of A CT for different values of X e (Fig. 3.14(a)) and in the opposite case (Fig. 3.14(b)). 

Summary of the results can be seen in Figures 3.15 and 3.16, where the surface of bi- 
furcation points is plotted against both the shape and energy biaxialities, in the rod-like and 
disc-like molecular regime, respectively. The line drawn on the upper surface (representing 
the transition from isotropic phase to Ny) in Fig. 3.15 is a family of Landau points where the 
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(a) Bifurcation temperature as function of shape 
biaxiality, for some A e . 




-0.2 -0.1 0.1 0.2 



(b) Bifurcation temperature as function of en- 
ergy biaxiality, for some X a . 



Figure 3.14: Concurring biaxialities and their influence on bifurcation diagrams for 
N[/ + - Nb transition for p* = 0.18. A CT varies along the thick, gradient line in 
Fig. 3.10, and A e is changed by increasing \e x — e y \ while keeping e z = 0.2. 



direct second order transition to biaxial nematic from isotropic liquid phase takes place. Those 
points will be taken into closer consideration in the following section. 



Landau points 

Previously, we have seen how the parameters A e and X a , representing the degree of overall 
biaxiality of the intermolecular potential energy (3.21), influence the bifurcation temperature 
and density. In present section we will locate the Landau points in (A CT , A e ) space, keeping in 
mind that X a is changed along the path drawn with gradient in Fig. 3.10, while in the case of 
X e we hold e z = 0.2 and starting from A e = 0.0 for e x = e y , vary the energy biaxiality by 
increasing \e x — e y \. We will also comment on the meaning of the hard biaxial ellipsoids self- 
dual plane (X s e d , X s a d ) (3.32)-(3.34), and its relation to Landau points in present model. Firstly, 
we will investigate uniaxial energy A e = 0.0 ((d) e in Table 3.3) and look for shape-induced 
effects, then we will turn to the case of fixed shape biaxiality and by changing A e pursuit the 
Landau point, and finally we will show how both A e and X a change the position of this point 
on the bifurcation diagram. 
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of shape (A CT ) and energy biaxiality (A e ). X a changes along the path in Fig. 3.10, A e is 
varied by increasing \e x — e y \ while e z = 0.2. Dark line on the upper surface marks 
the Landau points presented in detail in Fig. 3.18. 

In Fig. 3.17(a) we present the bifurcation line as function of shape biaxiality for \ e = 0.0 
for (d) e . It can be seen that there exists a point at the boundary of disks and rods, where the 
molecule belongs to neither one of those groups and Nu+ and Nu- phases are indistinguish- 
able. This is an isolated Landau point, located near X a = 0.7163 for (g) a (see Table 3.3). It 
differs from the predicted for hard biaxial ellipsoids [29, 80] self-dual geometry a x = sJ<J y a z 
(3.32), which in our parametrization would be at A CT = 0.6412 for (c) CT . As we can see, the 
attractive forces shift the position of Landau point towards higher values of \ a . 

If we can induce the transition directly to biaxial nematic phase by changing A^, then it 
should be possible to do the same thing for A £ . Figure 3.17(b) shows such bifurcation diagram 
as function of energy biaxiality for A CT = 0.70 for (/%. It is clear that for \ a ^ by changing 
A e we can induce a Landau point. We should recall, however, that we do not cross the bound- 
ary between "rod-like" and "disc-like" energies (in terms of (e x ,e y ,e z ) or A e ). It turns out 
that for shape biaxialities close to the ones giving Landau point by strengthening the lateral 
interactions we can also induce a self-dual point (e.g. the left one in Fig. 3.17(b)). 

Plot presented in Figure 3.17(b) gives us also some insight into the issue of dependence of 
Landau point position on A e and A^. There exist two such points in Fig. 3.17(b) for A CT = 0.70 
whereas in Fig. 3.17(a) for a bit higher shape biaxiality \ a = 0.7163 there is only one Landau 
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point visible. Given all this we can wonder how the location of this point changes under the 
varying biaxialities. We can plot Landau point position and divide the (A CT , A e ) space into disc 
and rod-like regions, as is shown in Figure 3.18, while keeping in mind the way (a x , a y , a z ) 
and (e x , e y , e z ) were changed and combined into X a and A £ . It can be seen that it is "easier" to 
procure a transition directly to biaxial nematic phase for the case of the same sign of X a and 
A e , whereas in the case of negative X e it takes stronger (relatively to A CT ) energy biaxialities to 
induce a Landau point. The increase of A CT resulted in meeting of the two Landau points visible 
on Fig. 3.17(b) at a single point for (A CT , A e ) = (0.7163, 0.0) as presented in Fig. 3.17(a). 

More interestingly, we can make some comments on the meaning of the self-dual condi- 
tion for hard, biaxial ellipsoids (3.32)-(3.34). We can see that it is possible to find Landau 
point away from the square root plane (X s /, X^f) (3.33)-(3.34), it is visible in Fig. 3.17(a). We 
have also investigated vicinity of the point lying in this plane at (X s e d , X^) = (0.175, 0.6412) 
((h) £ and (c) a in Table 3.3) and found the Landau point at (A e , A CT ) = (0.175, 0.6596) ((h) e 
and (d) a ). We can conclude, although based only on the four points (lying on the upper line in 
Fig. 3.18), that square root rule (3.32) is a good candidate for the starting position in the search 
for Landau points; we have shown that self-dual points can exist beyond (X s e d , X^) plane, but 
also in its vicinity. As we can see, the upper line in Fig. 3.18 crosses the vicinity of the point 




Figure 3.16: Bifurcation temperature t* for p* = 0.18 for Nj/_ - Ng transition as 
function of shape (X a ) and energy biaxiality (A e ) in the disc-like regime. X a changes 
along the path in Fig. 3.10, A e is varied by increasing \e x — eJ while e z = 0.2. 
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(a) A £ = 0.0 for (d) e (see Table 3.3). Landau 
point at \ a = 0.7163 for (g) a ((a x , a y , a z ) = 
(1.4, 0.544, 3.17)) is marked. Dotted line rep- 
resents the self-dual geometry for hard ellip- 
soids (3.32) from [29] at \ a = 0.6412 (c) CT 
({a x , Gy , a z ) = (1.4,0.637,3.077)). 



(b) X a = 0.70 (/) CT ; two Landau points are vis- 
ible at A e = -0.1974 and A e = 0.0662 for(see 
Table 3.3) (b) e and (/) e , respectively. They meet 
at Act = 0.7163 for (g) a at the Landau point in 
Fig. (a). 



Figure 3.17: Bifurcation temperature t* for p* = 0.18 as function of A CT (a), and X e 
(b). A CT varies along the path drawn in Fig. 3.10, X € is changed by increasing \e x — e y \ 
while e z = 0.2. 



on the plane in (A e , X a ) space defined by square root rule. When we consider the complicated, 
non-linear manner in which {e^} and {<Xj} contribute to the biaxial Gay-Berne potential, it is 
surprising that the Landau points can be found close to the plane obtained from the simple 
square root rule: e x = v /e^el, a x = ^fdyd~ z . 

The bifurcation analysis can give the position of Landau point (using Eq. (2.72)); however, 
in its vicinity it becomes increasingly hard to stabilize the uniaxial nematic reference phase 
because it looses stability at this point. Therefore, in order to confirm the places where self- 
dual point appears to exist, we have employed the method of minimisation of the Helmholtz 
free energy (3.5) with trial one-particle distribution function (3.9) by solving the equations 
(3.10) and obtained the temperature and density of the transition by calculating order parame- 
ters (2.49). An example of the behaviour of the leading order parameters at the Landau point 
is presented in Fig. 3.19, and for comparison we also show the plots away from this point in 
rod-like regime, both obtained from the minimisation method. 
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Figure 3.18: Division of (A e , A CT ) space by a line of Landau points separating disk- 
like and rod-like states. \ a runs along the thick, gradient path in Fig. 3.10, and \ e is 
changed by increasing \e x — e y \ while keeping e z = 0.2. The line joins the marked 
Landau points, and the dotted lines indicate the points where square root rule (3.32) is 
fulfilled for {e^} and {a;}, their intersection (point on the (X s /, \ s a d ) plane) is marked. 

3.4.5 Summary 

The uniaxial Gay-Berne potential [96, 97] achieved an amazing, compared to its simplic- 
ity, success in predicting liquid crystal phase behaviour. As pointed out before, most of the 
mesogenic molecules are not uniaxial. They are not biaxial either, but effectively in the uni- 
axial phase they can be approximated by ellipsoids of revolution, as was done in the original 
Gay-Berne interaction, and in biaxial phase by D 2 h - symmetric objects, as in the potential 
developed by Berardi, Fava and Zannoni [36]. We have investigated the latter model. 

We compared the bifurcation analysis of DFT in low-density approximation with the 
Monte Carlo results. As expected the bifurcation temperatures were significantly higher than 
those following from simulations. This tendency was more apparent in case for I so - ]%+ 
transition, while for N[/+ - N B bifurcation the approach provided more accurate estimates in 
relation to Monte Carlo. The comparison was only qualitative because the simulation study 
neglected long-range corrections, the inclusion of which makes the results of Monte Carlo 
closer to those of DFT [108]. 

Then, in order to provide some insight into the way the reference state influences the phase 
diagram, we calculated the bifurcation for Nu + - Nb transition using the minimisation of the 
Helmholtz free energy. This method does not involve the truncation of the expansion of the 
pair direct correlation, but uses the whole function, and therefore, compared to bifurcation for 
second order phase transition, provides the estimation of the influence of precisely calculated 
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Figure 3.19: Temperature behaviour of leading order parameters at fixed density p* = 
0.18 in the vicinity of Landau point (on the left) at \ a = 0.7163 ((g) & m Table 3.3) 
and, on the right, away from that point for rod-like molecule (giving Njj+) for A CT = 
0.58, X e = -0.06 ((a x ,ay,a z ) = (1.4, 0.714, 3.0), (e x , e y , e z ) = (1.7,1.0,0.2)). As 
obtained from the minimisation of the Helmholtz free energy along Eq. (3.5), (3.8), 
and (3.10). 



one-particle distribution function of reference phase. The results show that the inclusion of 
terms with higher angular momentum index than 2 in the expansion of pair direct correlation 
function has secondary meaning. 

Next, we continued to investigate the influence of introduced earlier shape biaxiality \ a 
and energy biaxiality A e on the bifurcation diagram, by defining a path in the six dimensional 
space of potential shape and energy parameters. We presented the diagrams for uniaxial cases 
of the potential (A e = 0.0) for different values of the biaxiality of shape, and uniaxial shape 
(A CT = 0.0) was studied for given energy biaxialities. We also addressed the issue of opposite 
signs of A CT and A 6 . With increasing positive A CT , for uniaxial - biaxial nematic bifurcation 
the line of transitions corresponding to A e < was lying in temperatures higher than the one 
associated with A e > 0, apart from the point when A CT became large enough to dominate the 
transition. The results showed how the Nu - Ng transition is influenced by the Gay-Berne 
interaction biaxiality originating from different sources. They also suggested that the model 
can exhibit a Landau point; this possibility was also investigated. 

We found that by increasing A CT for A e = 0.0 we can cross the boundary between the oblate 
and prolate nematic phases and locate the Landau point when the distinction between them is 
not possible, where I so, N;y„, N u+ , and N B phases "meet", and the second order isotropic - 
biaxial nematic phase transition occurs. This shape-induced self-dual point was found to be 
located at higher shape biaxiality than the one obtained for hard molecules. In terms of A CT , 
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the difference between the shape biaxiality at Landau point following from our analysis (with 
\ t = 0.0, so the biaxiality of the interaction originated only in shape of molecules), and the 
one obtained for hard biaxial ellipsoids was found to be 0.0751. This difference is the estimate 
of the degree of the influence of attractive, anisotropic forces on Landau point position. The 
self-dual points were also found to occur for fixed shape biaxiality when A e was changed. In 
this case we varied the energy biaxiality by going from the model of strong lateral interactions 
to the one where molecules are most attracted to their faces (we did not cross the "rod-like" 
and "disc-like potential" boundary in terms of A e ). When energy parameters fulfilled the self- 
dual square root rule for hard ellipsoids, the Landau point was found for the \ a different 
from the one predicted for hard potentials by 0.0184. The same difference, far away from 
£x = yj e y e z, for \ e = 0.0 was more than four times higher. The matter requires further stud- 
ies, but our results suggest that the meaning of the square root rule maintains its importance 
as the qualitatively accurate estimates for Landau region for soft biaxial Gay-Berne interac- 
tion. We concluded our findings by the division of (A e , A CT ) space into regions belonging to the 
oblate and prolate states. It proved that the line of Landau points crosses the vicinity of the 
point following from square root rule for {ej} and {o"i}. The existence of Landau points their 
position and order of the transition were confirmed by the minimisation of the Helmholtz free 
energy. 

The thermotropic biaxial nematic phase discovered in simulations by Berardi and Zannoni 
was reported for the model of stronger lateral interactions, i.e., where the side-to-side config- 
uration is preferred (see Fig. 3.5), that is, for opposite signs of shape and energy biaxiality. 
Our results suggest that while the temperature of N u+ ~ Ng transition is relatively higher for 
negative X e (while \ a is greater than zero, see Fig. 3. 13) it is necessary to provide significantly 
lower (negative) energy biaxialities to induce the Landau point in that case (see Fig. 3.18). 
Whereas for positive biaxialities the transition directly to biaxial nematic can occur for mod- 
erate values of A e and A CT . We shed some light on the issue of the case A CT > 0, A e < in the 
last chapter, where we take into account uniaxial and biaxial smectic-A phases; it is natural to 
suspect that those structures are preferred over N B when the interaction favours face-to-face 
configuration over side-to-side. 

Apparently, in case of biaxial nematic, the significance of the Gay-Berne interaction is 
much smaller than for uniaxial phases. Therefore, in order to continue the search for factors 
behind the stabilisation of N B , we need to study a model that can be more closely related to 
the real systems. We turn our attention to the bent-core molecules. 
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3.5 Systems of bent-core molecules 



In present section we consider the models of bent-core (banana) molecules built from 
Gay-Berne segments 3 . We study shape-related effects and the influence of dipole-dipole in- 
teraction. Firstly, we take into account the case of D^h - symmetric arms and join two and 
three prolate, soft ellipsoids, as presented in Fig. 3.20, to model banana- like molecules with 
bend (opening) angle 7. The system is studied for fixed density, and the resulting bifurca- 
tion diagrams are presented in the plane of temperature t* (see (3.20)) and 7. We especially 
focus on Landau point position, as obtained from (2.72), and in order to shed more light on 
this issue, we also take a look at the behaviour of Landau region when the arms deviate from 
Dvoh symmetry; we incorporate two soft, biaxial ellipsoids into a bent-core molecule using 
the Berardi-Fava-Zannoni potential [36], studied in the previous section. Inspired by recent 
developments [46], we also comment on the asymmetric bent-core molecules with arms mod- 
elled by soft, Dooh - symmetric ellipsoids of different elongations, using potential from [103]. 
Finally, we turn to the effects of polarity, and we include dipole-dipole interaction between the 
dipole moments located along C 2 symmetry axis of the molecule (see Fig. 3.20). We present 
the phase diagrams in (£*, 7) plane for fixed density at bifurcation point, again focusing on the 
behaviour of Landau point. 

We begin with a description of the model, then we present the bifurcation diagrams for 
non-polar case for uniaxial and biaxial arms, and then we study the influence of dipole-dipole 
interaction. 



3.5.1 Models of a bent-core molecule 



This section is devoted to the description of the pair potential that was used as an inter- 
action between two bent-core molecules. We are using two already mentioned models, those 
are: the uniaxial Gay-Berne [96, 97] and its biaxial version as proposed by Berardi, Fava, and 
Zannoni [35]. Fig. 3.20 shows the construction. Each of the parts interacts with every one on 
the other banana, while the interaction between segments within single molecule, giving an 
irrelevant constant, is disregarded. We consider rigid bananas, i.e., we do not allow the change 
of the angle 7 or dimensions of the ellipsoids during the calculations. 



Which is a simple way to do it, but at least guarantees the basic correct behaviour of the molecule. 
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The intermolecular potential can be written as (see Fig. 3.20): 




ab) i 



(3.35) 



where, as before, Vqb is either the uniaxial Gay-Berne potential, described in Sec. 3.4.1, or its 
biaxial version studied in previous section and introduced in Sec. 3.4.2, Ujj, [i = 1 . . . 2, j = 
1 ... 3) are unit vectors parallel to the Gay-Berne segments symmetry axes, are vectors 
connecting the centres of molecules. For D 2 h - symmetric Gay-Berne, in order to describe 
the orientation of the ellipsoid, it is not enough to give the vector parallel to the longer axis, in 
this case we fixed the remaining axes with respect to the plane of the angle 7, namely one axis 
was set to lie in this plane and the other to be perpendicular to it. For the uniaxial Gay-Berne 
the empirical exponents v, \x are chosen to be to 1 and 2, respectively, and the elongation of 
the ellipsoids is 5 : 1, and anisotropy of the potential 4 : 1, in terms of constants from (3.16) 
and Eq. (3.19) it means k = 5 and x = 1/3 (see Fig. 3.4(b)). The molecular elongations are 
chosen to correspond to the dimensions of the bent-core compounds. 

We have also added the dipole moment localized in the fixed position, as shown in Fig. 3.20, 
along molecular C 2 symmetry axis. The dipole-dipole interaction used is of a well known 
form: 



Vdd (a2i,/*2> ? 



& • /2 2 - 3 (/It • f ) (gg_j) 




Figure 3.20: Construction of bent-core molecule and pair potential. On the left there 
are frames of reference associated with banana 1: h\, b^, 63, and 2: b^, b|, b| (h\ 
and bj are perpendicular to the picture surface). 
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Table 3.4: Dipole-dipole interaction contribution to total potential energy in the 
ground state. 

/i* \V DD /{V b3 + V DD )\ 



2.8 


0.50 


2.2 


0.40 


1.6 


0.25 


1.5 


0.22 


1.2 


0.15 



where \i i = an d where fi i is a unit vector, and jj* is the dimensionless (see (3.20)) 

magnitude of the dipole moment. To provide some way of relating ji* to the strength of Vqb, 

in the ground state for given values of fx* (see Table 3.4). This 



VPD 



we have calculated 

ratio gives relative weight of contribution from dipole-dipole interaction to total potential; it 
will be used as a measure of dipole strength. 

We employ the units of energy (e ) and distance (a ) as described by (3.20). Results pre- 
sented in this section follow from using in place of V in Eq. (3.1)-(3.2) Vm or V&3, with Vdd 
added in polar case. 

The absolute minimum of the potential between two banana-like molecules (the ground 
state) for the non-polar uniaxial model corresponds to the configuration where h} and bf 
axes of frames associated with molecule 1 and 2 (see Fig. 3.20) are parallel and the vector 
r is oriented along b\. The inclusion of the dipole-dipole interaction changes that picture; 
in the ground state relative orientation is the same, only r is parallel to . In numerical 
analysis for two-part banana, to prevent the infinite values of potential when dipoles are lo- 
calized at the same point, we have introduced a hard sphere between the arms with radius of 
0.5 do. Fig. 3.21 (not containing the just-mentioned sphere) shows the equipotential surfaces 
for the model where uniaxial Gay-Berne potential was used, for four relative orientations of 
two molecules. 

Every bifurcation diagram presented in this chapter was calculated at fixed density chosen 
in the following way. If we calculate the inverse of the excluded volume for parallel bent-core 
molecules we will receive values of 0.024 a^ 3 and 0.016 <7q 3 for two, and three -arm bananas, 
respectively. The densities on the diagrams were set to be of that order, which corresponds to 
packing fractions (molecular volume per average volume of molecule) of 0.2 and 0.3, being 
slightly lower than typical values in uniaxial nematic which range from 0.4 to 0.6. 

In the following section we present the bifurcation diagrams for bent-core molecules with- 
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Figure 3.21: Sample equipotential surfaces for banana composed of uniaxial Gay- 
Berne molecules, 3 parts (on the top) and 2 parts (on the bottom), for opening angle 
of 0.7n = 126°. On the right, case with dipole-dipole interaction with dipole strength 
H* = 2.0. On the left, case without dipoles. Surfaces show the potential equal and 
-0.2. 

out dipole moments. 

3.5.2 Shape-induced effects 

In present section we seek the answer for the most straight forward question of what hap- 
pens to the bifurcation temperature as we change the opening angle 7, and also whether we 
can find Landau points, where direct transition from the isotropic to biaxial nematic phase 
takes place. Since we pursuit only the shape related effects, we do not introduce any dipole- 
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(a) Diagrams for two-part, non-polar banana, for two values of bifurcation density p* 
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(b) Diagrams for non-polar model of three parts, for two bifurcation densities p*. 

Figure 3.22: Bifurcation diagrams for model bent-core molecules without dipole- 
dipole interaction. The temperatures are scaled by tg , the bifurcation temperature for 
7 = 90° (see Table 3.5). 



dipole interactions yet. Firstly, we present the results for model where the arms of bent-core 
molecule are uniaxial ellipsoids of revolution, and then we take a brief look at the case of D 2 h 
- symmetric building blocks. 
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Uniaxial arms 

We can think of the arms of a banana (see Fig. 1.13) as of uniaxial ellipsoids. This ap- 
proximation probably will not reflect all properties of the bent-core molecule correctly, but 
can serve as a starting point for more complex models. 

We have taken into account the molecules composed of two and three Gay-Berne inter- 
acting parts. In Fig. 3.22 we show the resulting bifurcation diagrams for non-polar banana 
for two densities. As can be seen, the opening angle at Landau point, for which the direct 
transition from isotropic phase to the biaxial nematic takes place, for two-part banana is equal 
to 107°, which is in agreement with the results for hard interactions [38] and close to the mean 
field predictions (109°) [39]. Interestingly, for bent-core molecules composed of two parts 
the correction to the position of the self-dual point from attractive forces is vanishingly small. 
We remember that it was not the case for convex biaxial molecules (see Fig. 3.17(a)). The 
model of three arms, as is shown on Fig. 3.22(b), has the Landau point localized near the right 
angle, which can be considered to be in some agreement with observations in [61], where the 
mesogenic substance consisting of molecules with bend angle of 90° was found to give rise to 
biaxial nematic phase. The summary of Landau points positions is presented in Table 3.6. 

In order to pursuit the behaviour of Landau point in the two-arm model, we have also 
taken into account asymmetrical bananas, i.e., we constructed the bent-core molecules from 
two uniaxial ellipsoids with different elongations. We used the extended Gay-Berne potential 
[103] and performed the calculations for a range of asymmetries, up to a point where one arm 
was four times longer than the other. No significant differences in the position of the Landau 
point were found, it was located, as for symmetric bananas, at 107°. 

Non-polar model of uniaxial arms is able to produce an isolated self-dual point, given in 
Table 3.6. Following sections show two models where that behaviour is different, namely, 
there appears a line of Landau points. 

Table 3.5: Temperatures of bifurcation from isotropic phase for 7 = 90°. 



Model 


P* 


t* 




0.04 


0.84 


two-arms bent-core molecule 




0.06 


0.96 


three-arms bent-core molecule 


0.026 
0.04 


1.34 
1.71 
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Table 3.6: Landau point position versus strength of the dipole p*. 



Bent-core model 




opening angle at Landau point 




0.0 


107° 


two-arms 


1.2 


104° 




1.5 


103° 






p* = 0.026 p* = 0.04 




0.0 


89° 89° 




1.2 


86° 86° 


three-arms 








1.6 


74° - 86° 74° - 86° 




2.0 


63° - 86° 63° - 80° 




2.8 


83° - 97° 82° - 92° 


two biaxial arms 


0.0 


121° - 128° 



Bent-core systems composed of biaxial parts 

It is of interest to take into account the model where the arms of a banana are not uniaxial. 
One possibility would be to consider the bent-core molecules built of the biaxial parts. We 
have done so, and present the resulting phase diagram in current section. 

We take the biaxial Gay-Berne ellipsoids as building blocks. They are oriented in the 
following way. We fix the plane of the longer axis, to lie in the plane of the bend angle, 
while the shortest axis is chosen to be perpendicular to it. In this way we construct the bent- 
core molecule as shown in Fig. 3.20, and build the pair interaction by substituting V GB in 
(3.35) with the biaxial Gay-Berne potential (3.21), defined in Sec. 3.4.2. The resulting pair 
interaction has the global minimum in the same configuration as the uniaxial two-arm model 
described in the previous subsection. 

We are mainly interested in the position of the Landau point in this model. We have 
chosen the shape of constituting biaxial ellipsoids to be (a x , a y , a z ) = (1.2, 0.514, 3.4) and 
the potential parameters (e x , e y , e z ) = (1.0, 1.4, 0.2). As mentioned in the previous chapter 
those parameters make the attractive forces strongest in the face-to-face configuration, that is, 
in the direction of the shortest axis (see Sec. 3.4.2). This model exhibits molecular biaxiality 
in the limit of 7 = n and 7 = 0; however, it is as well clear that there should exist an angle for 
which Landau point appears. Indeed, the numerical analysis shows the existence of self-dual 
region, in the form of the line of Landau points, as can be seen from Fig. 3.23. For the biaxial 
arms interacting by biaxial Gay-Berne potential a line of I so - Nb transitions is possible. It 
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Figure 3.23: Bifurcation diagram for bent-core molecule constructed from biaxial 
ellipsoids interacting by Berardi-Fava-Zannoni potential [35, 36], for fixed density 
p*. Landau line is in range of opening angle 121° < 7 < 128°. 

starts near 7 = 121° and ends for 7 = 128°, and reduces to a single point with a decrease in 
arms biaxiality. 

The matter should probably be investigated further, since the Landau region obtained for 
the model with biaxial arms is closer to the experimental findings, not to mention that the 
topology of the diagram is qualitatively different from what was published in the literature so 
far. 



3.5.3 Polar case 

As was mentioned before, molecules studied in [57, 58, 61], i.e., those for which the bi- 
axial nematic was observed, are polar. Therefore, it is interesting to look at the effects of 
dipole-dipole interaction on stability of the biaxial nematic phase in model bent-core systems. 
We have introduced the dipole moment along the molecular C2 symmetry axis, as depicted in 
Fig. 3.20. In current section we present the bifurcation diagrams for two and three uniaxial 
arm models for a couple of dipole strengths. We again focus on the dependence of Landau 
point position on the dipole magnitude (see Table 3.6). 

Fig. 3.24 shows the bifurcation diagrams. As can be seen (Table 3.6) the inclusion of the 
dipole-dipole interaction shifts the Landau point for the model of two uniaxial arms towards 
lower angles (Fig. 3.24(a)). In case of bent-core molecule modelled by three uniaxial Gay- 
Berne parts (Fig. 3.24(b)), for values of dipole strength less than 20% of the total potential in 
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(a) Diagrams for two-arm model for two densities p* and three dipole magnitudes p*. 




(b) Diagrams for three-arm banana, for two densities p* and five dipole strengths p*. 

Figure 3.24: Bifurcation diagrams for systems of polar bent-core molecules com- 
posed of uniaxial, soft ellipsoids. Temperatures are scaled by tg , the bifurcation 
temperature for 7 = 90° in non-polar case (see Table 3.5). 



the ground state (see Table 3.4), the same tendency is maintained, i.e., Landau point is shifted 
towards lower angles up to a value of 86°. However, for dipoles above fi* = 1.4 (20% of 
total potential) isolated Landau point changes into the second order I so - Ng phase transition 
line. The range of this line measured in opening angles widens; for // = 1.6 it is 12° and 
for fx* = 2.0 exceeds 20°. For moderate density, its upper limit does not change and is at 
7 = 86° as long as the dipole strength /1* < 2.1 (38% of total energy). Then, the Landau 
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Figure 3.25: Landau regions as function of dipole-dipole interaction strength contri- 
bution to total potential (see Table 3.4). 

region begins to shrink and is shifted towards higher angles. The highest dipole studied is the 
one of n* = 2.8 for which Vdd is 50% of total potential. As can be seen from Fig. 3.24(b), 
the line of I so -N B transitions in that case is still getting shorter, and moves in the direction 
of higher angles. Fig. 3.25 shows the evolution of Landau region for varying dipole-dipole 
potential contribution to total energy. When Vdd contribution is more than 35% and less than 
38%, the Landau line attains its maximal width in opening angle. 

We have presented the diagrams for two densities selected from lower region, as described 
above. It can be seen (Table 3.6) that for strongest dipoles studied (p* > 2.0) some differences 
appear with varying density. Namely, the line of direct isotropic - biaxial nematic transitions 
gets shorter for higher density. 

Our results suggest that there exists a certain range of dipole strength which makes the 
appearance of biaxial nematic phase most probable. Based on the diagrams presented in 
Fig. 3.24, and the assumed model parameters, we estimated this region to be between fi* — 1.4 
and 2.1, that is, when dipole-dipole interaction contributes in range of 20 and 38% to total po- 
tential energy. 

Finally, we should mention that the only other studies in this direction are Monte Carlo 
simulations of lattice Lebwohl-Lasher model presented in [48], where the introduction of 
dipoles resulted in the line of direct isotropic - biaxial nematic transitions. 

3.5.4 Conclusions 

We studied the bent-core molecule model. Firstly, we took into account the bananas com- 
posed of two and three uniaxial Gay-Berne ellipsoids. We discussed briefly the case where 
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arms of the molecules interact via biaxial Gay-Berne potential, and also where one arm is a 
longer uniaxial ellipsoid than the other. The dipole-dipole interactions were introduced next, 
and the influence of the dipole strength on position of Landau point was studied. Our results 
qualitatively support the previous findings [38, 39], and show similar behaviour as the less 
realistic polar, lattice model [48]. 

The bent-core Gay-Berne model of two uniaxial arms exhibits the Landau point in agree- 
ment with previous studies for hard bananas in Onsager approximation and mean field [38, 39]. 
Its position was not changed significantly by the introduction of attractive forces, nor by the 
inclusion of anisotropy of the molecule. In case of the three -Dooh - symmetric arms, opening 
angle near the right angle was distinguished as the one where Landau point appears. Recent 
experimental studies [61] are in line with these findings. 

We tried to address the problem of disagreement between the experimental estimates of 
the opening angle in the biaxial nematic phase (predicted to be near 140°) and the theoretical 
results by switching to biaxial arms interacting via biaxial Gay-Berne potential. In that model 
with the increase of arms biaxiality, Landau point evolved to the line of direct I so - Nb second 
order phase transitions in range of opening angle 7 between 121° and 128°. Probably, inclu- 
sion of dipoles or increase in the arms biaxiality or taking into account molecular flexibility 
could widen that region and possibly shift it towards higher angles. 

That behaviour was not possible to acquire in non-polar case of uniaxial arms. The line of 
direct I so - Nb transitions appeared in the model of three-part banana with dipole localized 
along the C 2 symmetry axis. It appeared that it is possible to provide dipoles strong enough to 
generate the line of Landau points, which widened with the increase of the dipole strength up 
to its maximal range in the opening angle (63° < 7 < 86°), when the dipole-dipole interac- 
tion constituted between 35 and 38% of total potential. Stronger dipoles began to shrink that 
line and move it towards higher angles. That should provide some hints for future research. 
Possibly a study on the orientation of the dipole as well as on the shape of arms and shape of 
potential is in order, for instance, we could wonder what would happen when shape of arms is 
close to the self-dual geometry (found in Sec. 3.4.4). 

Our studies clearly indicate that strong dipoles were a necessary feature for experimental 
discovery of the biaxial nematic phase in bent-core systems. In polar case we saw that the 
dipole-dipole interactions not only were able to shift the position of the Landau point (move 
it towards lower angles) but also change the topology of the bifurcation diagram to include 
the self-dual line of direct isotropic - Nb transitions. We can speculate that if the similar 
mechanism is at work in the substances where the biaxial nematic was observed, then the 
dipoles indeed became crucial, not only increasing the biaxiality of a molecule but also shifting 
the N B region to lower angles where the realization and observation of the elusive phase 
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becomes possible. 

We were unable to account for the experimentally predicted value of the opening angle 
of 140°, as observed in the first class of the bent-core biaxial materials. The stability of 
the biaxial nematic for the angles of this magnitude, in comparison to Iso and phases, 
was very low. The introduction of dipole-dipole interaction shifted the Landau point in the 
opposite direction, so probably it is safe to conclude that the polarity alone cannot explain the 
disagreement between theoretical and experimental results. Only the inclusion of biaxial arms 
interacting via biaxial Gay-Berne potential gave the values of opening angle somewhat closer 
to 140°. 



Chapter 4 
Preliminary results 



The knowledge of phase diagram does not give the complete information about behaviour of 
a system, e.g., from the knowledge of the density-temperature region of stability of a given 
structure we cannot extract the information about response of a system to a given deformation 
of the director field. For many applications this kind of information is crucial. This chapter 
is devoted to the two issues that were disregarded in the preceding part of this thesis. One 
of them is the investigation of a phenomena related to small distortions of the homeotropic 
biaxial nematic phase and description of such state in terms of the biaxial elastic constants. 
The other is the matter of spatially non-uniform liquid crystal states. Presently we show the 
preliminary results concerning those issues. 



4.1 Biaxial elastic constants 



Response of the biaxial nematic state to a small perturbations is usually described by the 
contribution to the free energy coming from the distortion of the directors, expanded in terms 
of their gradients. The coefficients of this expansion are called elastic constants, and are of 
great importance to both theoretical and experimental studies. They influence most of the 
observed phenomena in liquid crystals. To name only a few, we can mention their role in 
light scattering effects, response to external fields [64], shape of disclination defects [109], 
and nematodynamic flow [110]. The elastic constants are crucial quantities in applications 
of liquid crystals in displays [111], and also are of importance in polymer-dispersed systems 
[112]. They are as well significant for the stability of the fhermotropic biaxial nematic phase. 
The low relative magnitude of elastic constants corresponding to the biaxial distortions in 
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comparison to uniaxial ones can be one of the reasons for difficulties in stabilizing N#. 

The usual way of approaching the issues of elasticity is to study the elastic free energy 
density dependence on the local deformations. The first derivation is due to Oseen [113], 
later phenomenologically acquired by Frank [109], and also by Zocher [114]. In those classic 
papers the liquid is studied by means of the director field n(r) and the deformations of this 
field give rise to the increase of free energy 1 , typically described by expansion in gradients 
of n. In the case of biaxial nematic phase, we have two additional directors denoted by m 
and 1 which are not independent on each other and on the uniaxial director. The fact that the 
directors form an orthogonal tripod means that we cannot perform an independent distortion 
of one director field without affecting other two. That, in turn, leads to the impossibility of the 
expression of the elastic free energy density as a sum of independent deformations, as in the 
uniaxial case. 

Current section by no means is a complete study on elasticity issues; it serves merely as 
an introduction to further research. We are using the biaxial Gay-Berne potential, studied 
in the previous chapter, and present the set of elastic constants in the case of biaxial prolate 
molecular shape, and in the vicinity of Landau point. We follow and employ the approach 
from [115]. We should note that there are other known distortions of the director field where it 
looses the continuity, these are called topological defects [116, 117] and we are disregarding 
them altogether. 



4.1.1 Elastic constants in the biaxial Gay-Berne model 

The results presented in this section follow from the application of the formulas from 
[115]. We will not repeat the full derivation, since it can be found in the literature. Instead we 
briefly write down the basic formulas. 

The elastic free energy density for biaxial nematic, in absence of polar and chiral ordering, 
can be written in many ways. It is not possible to derive a formula where the coefficients of the 
expansion in gradients of n, 1, and m are independent, due to the fact that the three directors 



At the level of the molecular interpretation the free energy density of the small elastic deformations is always 
positive. 
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form the orthogonal tripod. We use the form developed in [1 18]: 

^elastic =^i(divl) 2 + K l2 (\ ■ rotl) 2 + K m (\ x roti) 2 + 

+K ml (div m) 2 + K m2 (m ■ rot m) 2 + K m3 (m x rot m) 2 + 

+K nl (div fi) 2 + K n2 (n • rot n) 2 + K n3 (n x rot n) 2 + 

+K mn (m ■ rot fi) 2 + K n i(h ■ rotl) 2 + K lm (\ ■ rot m)+ ( J) 

+Ki A 6\\ (1 ■ V)l — 1 div 1 + A' m4 div (m ■ V)m — m div m + 

+i^ n4 div (fi ■ V)n — fi div n , 

where Kn, K i2 , K i3 for i = n,l,m denote the elastic constants corresponding to the deforma- 
tions of director field fi, 1, rh called splay, twist, and bend, respectively. K mn , K ni , and Ki m 
denote the constants associated with mixed distortions, and the last three terms are propor- 
tional to full divergences, and therefore play a role only close to intrinsic surfaces (defects). 
Thus analysing the bulk effects we are concerned only with 12 elastic constants. They are not 
independent, and some of them may be negative, but the overall F e i ast i C > 0. There are other 
methods of expressing the elastic free energy density, see e.g., [119]. 

In order to obtain the molecular expressions of elastic constants, we recall the expansion 
(2.26) of the free energy around a given reference state g re f now representing the undeformed 
structure. In order to evaluate (2.26) so it contains the local terms, which will allow to re- 
late the coefficients Kij in (4.1) to microscopic quantities, we perform the following gradient 
expansion of the one-particle distribution function of the deformed state £>(r, CI) [115]: 

f?(r 2 , Q 2 ) = £?(?!, Q 2 ) + r« 2 d a £?(r l7 n 2 ) + ^r&rfAfy n 2 ) + . . . , (4.2) 

where Fi 2 stands for separation vector, r" 2 denote its Cartesian components, and summation 
goes over the repeated indices. Now we can insert the above to (2.26) and obtain [115]: 

(3Felastic = J d Z Y f3¥ e l as tic = J ^%|;" Q/ 3 7 ' 7 " (^y 'j^Tiy, ') , (4.3) 

where we have set ^(r, ft) = g(Cl^ ■ h^ k \r)) with n^'^ denoting the kth director of the dis- 
torted set, i.e., we have assumed deformed director fields to be slowly varying (see [120]), and 
again the summation convention is used. In the above 

M k'k" a ^'j" = \J dr^dVtxd^c^tl^tl^rn, [e])rf 2 rf 2 A fc /y (OOA^y,^) , (4.4) 
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where 

A w (fi) = -2^V (4.5) 
diljk 

Now by comparing (4.1) and (4.3) we can express Ky in terms of M k > k » a/3 > » . For exact 
derivation please refer to [115]. The results presented in current chapter were obtained by 
calculating the integrals in Eq. (4.4) in low-density approximation (3.1), i.e., by setting 

c 2 (n u fi 2 , r 12 , [q]) = exp[-/3 U{Q U fh, ?i 2 )] - 1 , (4.6) 

where U is defined in Sec. 3.4.2, Eq. (3.21), and then using Eqs. (B1)-(B15) from [115] to 
obtain the elastic constants. We have also used the expansion of the one-particle distribu- 
tion function in order parameters (2.48) remembering that g(f, f2) = pP(Cl). To calculate 
the derivatives of one-particle distribution function in (4.5) we employed Eqs. (A1)-(A4) and 
Eq. (A9) from [115]. 

In [115] set of elastic constants was calculated for a lattice model of biaxial molecules 
[85] interacting via a simplified version of two-point potential (3.11) (with t> ,o = ~1> v o,2 = 
v 2t o = 2 A, and v 2j2 = —2 A 2 , which is a version of the dispersion model studied in [86]). The 
splay-bend degeneracy was proven and observed for all directors. The authors also found one 
of the biaxial elastic constants to be negative for every parametrization studied, but with the 
smallest absolute value 2 . The relative values of the elastic constants associated with uniaxial 
director were found to be greater than biaxial ones. The two became relatively closer at the 
Landau point, which is believed to exists at A = 1/ \/6 in that model [88], but still K n \, K n2 , 
and K n3 where higher than their biaxial counterparts. 

We present the plots of the set of elastic constants for a biaxial case calculated within low- 
density approximation for the biaxial Gay-Berne model. As the results here are considered 
preliminary and a more detailed study is planned, we focus on the 9 bulk elastic constants 
Kn, K i2 , K i3 , for i = n,l,m. For comparison with the previous studies and to see how the 
biaxial constants behave close to second order I so - N# transition at Landau point, we choose 
to study the geometry where we found the Landau point (visible in Fig. 3.17(a)). Also to shed 
more light on the case of opposite sign of energy and shape biaxiality, we study two cases for 
rod-like biaxial molecule with A e > and \ e < 0. We carried out the calculations for three 
sets of shape and energy biaxialities (3.31). Firstly, we take into account the shape studied in 
Monte Carlo simulations [35]: (a x , a y , a z ) = (1.4, 0.714, 3.0) (A CT = 0.58) for two values of 



2 That is: the negative contribution to the free energy density from the deformation of director 1 was overcame 
by the effects of remaining, positive elastic constants associated with m and n. 
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X e = 0.04 ?L E = -0.04 X £ = 000 




1.6 1.8 2.0 2.2 2.4 1 6 1 8 2 .0 2.2 2.4 3.6 4.0 4.4 4.8 

t* t* t* 

Figure 4.1: Set of elastic constants as function of temperature, for biaxial Gay- 
Berne interaction. On the left the plot for model of positive A e for (e x ,e y ,e z ) = 
(1.0, 1.4, 0.2) (A e = 0.04) in the middle the negative case X e = -0.04 for (e x , e y , e z ) = 
(1.4, 1.0, 0.2), both for the shape studied by Berardi and Zannoni in [35] for A CT = 
0.58. On the right behaviour of the elastic constants in the vicinity of the Landau 
point for X € = 0.0 and A CT = 0.7163 (as seen on Fig. 3.17(a)). 



\ t = -0.04 ((e x , e y , e z ) = (1.4, 1.0, 0.2)) and A e = 0.04 ((e x , e y , e z ) = (1.0, 1.4, 0.2)) (the dia- 
gram for that model is presented in Fig. 3.8). Then we turned to the self-dual shape discovered 
earlier (see Fig. 3.17(a)) where (A CT , A e ) = (0.7163, 0.0) for (a x , a y , a z ) = (1.4, 0.544, 3.170) 
and (e x ,e y ,e z ) = (1.2,1.2,0.2). The calculations were performed for fixed dimensionless 
density p* = 0.18. 

The results are presented in two Figures: 4.1 and 4.2, where the reduced elastic constants 
K*- = (3 Kijp^e^o-Q 6 axe plotted against dimensionless temperature t*, where the units are 
set according to (3.20). The first one shows all the elastic constants calculated, and the second 
one presents a comparison of the coefficients associated with the deformation of biaxial direc- 
tors for A e > and A e < 0. 

The well known fact of equality of I<n and K i3 for i = m,nj (called a splay-bend de- 
generacy for uniaxial director) was found to be in force for all of the directors, which follows 
from the lack of director dependence in low-density approximation for pair direct correlation 
function C2 [121], and is due to keeping only leading order parameters in formula (2.48). It 
can be immediately seen from Fig. 4.1 that one of the elastic constants is always negative, this 
is a reflection of the impossibility of performing the distortion of only one director from the 
orthogonal set without affecting the other. The values of the constants corresponding to the 
biaxial directors are relatively low as compared with the uniaxial ones, both for opposite and 
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same sign of the energy and shape biaxiality. Close to the uniaxial - biaxial nematic transi- 
tion the biaxial constants became equal, as can be seen in Fig. 4.2. The situation changes 
qualitatively when we approach the Landau point. There the splay-bend degeneracy evolves 
to the equality of all the constants for biaxial directors 1 and m: Ku = Ki 2 = and 
K m i = K m2 = K m3 , and the latter become dominant. 

Fig. 4.2 shows the biaxial constants for two cases of the biaxiality of energy: \ e = —0.04 
and \ t = 0.04. It can be seen that although for \ e < the constants are higher there is no new 
quality revealed for the case of opposite signs of \ a and A e . 



4.1.2 Summary 

We have calculated the set of elastic constants for the biaxial Gay-Berne model in the low- 
density approximation. Taken into account the importance of the elastic issues, it seemed only 
fair to present the above study to give some insight into the effects of elastic deformations 
in the biaxial nematic phase formed by the molecules interacting via the biaxial Gay-Berne 
potential [36]. 

As we can see from Fig. 4.2, the models of opposite and same sign of shape and energy 
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Figure 4.2: Comparison of biaxial elastic constants for positive and negative energy 
biaxialities. 
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biaxiality do not give significantly different elastic constants. The previous analysis of bifur- 
cation diagrams was also unable to account for the apparent diversity of those cases, which 
was discovered in Monte Carlo simulations [35], its sources, therefore, must lie elsewhere; we 
address this issue in the next section. 

It is possible that the constants associated with biaxial directors are, comparing to uniaxial 
ones, too small and thermal fluctuations cause the lack of biaxial nematic ordering in experi- 
ment. The above results show that in the vicinity of Landau point at (A CT , A e ) = (0.7163, 0.0) 
the ratios of the elastic constants significantly change. For self-dual X a , at temperature 10% 
lower than the t* of Iso - N B transition they are K n : K ml : K nl = -0.71 : 2.46 : 1.0 (at 
t* = 4.14) as compared to K a : K ml : K nl = -0.18 : 0.33 : 1.0 (at t* = 1.65; 10% lower 
than t* for N u+ - N B transition) in rod-like regime for \ a = 0.58 and \ e = —0.04. As we can 
see one of the biaxial constants associated with biaxial director is dominant at self-dual point. 
However, a more detailed study, possibly involving bent-core molecules is needed. 

In the pursuit of the factors playing a role in stabilisation of the biaxial nematic, we con- 
tinue to the next section, where we make some comments on the simplest smectic phases. 

4.2 Orthogonal smectic phases 

As we have seen, both theoretical and experimental studies predict a stable biaxial nematic 
for bent-core molecules. However, in those systems the phase diagrams are dominated by 
smectic structures. Also one possible threat for the experimental observation of the spatially 
uniform biaxial phase comes from the strong possibility that the spatially non-uniform liq- 
uid crystal state can gain stability instead. Therefore it is important to study the competition 
between smectic and nematic phases. It seems, therefore, important to include those in the 
analysis. This section is intended as an introduction to a more complete study since the equa- 
tions presented here need to be solved and analysis of whether the smectic phases win over 
biaxial nematic needs to be conduced in more detail. 

There is a great variety of known smectic structures, and accordingly great interest is 
given to the investigation of materials that exhibit this kind of behaviour. Without trying to be 
exhaustive, we can mention the problem of the tricritical point on the line of uniaxial nematic- 
smectic-A transition, existence of which is still a matter of discussion, after over 32 years of 
research (see e.g. [77, 122]), and a rich number of higher ordered smectics, known as B n 
phases (n — 1, ... ,7), discovered in bent-core systems (see, e.g. [67]). 

In present section we briefly show the extension of the previous formalism to the case of 
the simplest spatially non-uniform liquid crystal states: the orthogonal uniaxial and biaxial 
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Figure 4.3: Picture of uniaxial smectic-A phase. The direction of density modulation 
z is parallel to the nematic director n, and perpendicular to the layers of thickness d. 



smectics. We start with an example of the order parameters calculated in L = 2 model includ- 
ing the smectic phases. Keeping in mind that the biaxial nematic phase range can be severely 
limited by the spatially non-uniform states, we again address the issue of opposite signs of 
energy and shape biaxiality (stronger molecular lateral attractions). As we mentioned, only 
for this model the simulations [35] discovered stable N B - We show how the range of biaxial 
nematic phase depends on the relative sign of \ e and \ a (defined in (3.31)). Then, as an intro- 
duction to a more complete study, we show the additional bifurcation formulas for 5 transitions 
involving the layered phases, these are nematic biaxial - smectic biaxial, nematic uniaxial - 
smectic uniaxial, nematic uniaxial - smectic biaxial, isotropic - smectic biaxial, and isotropic 
- smectic uniaxial. The numerical evaluation of these equations is yet to come. Firstly, we 
describe the basic smectic phases, and then show the usual methods used in the description of 
those structures. 

The simplest smectic phase is the so-called smectic-A (5mA). It is characterized by a 
long-range ordering of the orientational degrees of freedom, in the same way the nematic 
phase is; however, the positions of molecular centres of masses are not distributed randomly; 
there exists a tendency for them to on average align in layers perpendicular to the nematic di- 
rector n, while in each layer there is no long range translational order. Therefore a modulation 
of the density is visible. The inversion symmetry n — > — n and invariance under the rotation 
around n are upheld. Snapshot of such structure is pictured in Fig. 4.3. Since the layers are 
perpendicular to the director, SmA belongs to the class of orthogonal smectics; and due to 
the rotational symmetry around n this state is considered to be of uniaxial symmetry, and will 
be referred to by SmAu. It is possible that the biaxial nematic phase acquires the tendency to 
form orthogonal smectic layers. In that case the state becomes smectic with three directors, 
i.e., a biaxial smectic-A (SmA B ). Interestingly, the first historically discovered biaxially or- 
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dered phase was smectic [123]. Presently, we will only consider biaxial and uniaxial 5mA, 
for a description of layered structures with higher degree of order see, e.g., [64]. 

In the smectic regime the one-particle distribution function cannot depend only on orienta- 
tional degrees of freedom. To account for the average tendency of centres of mass to align in 
layers, it has to take, at least, one spatial direction as an argument. Traditionally it is denoted 
by z, and because it is chosen to be parallel to the uniaxial director n: z = zh we can conclude 
that P(Q, z) = P(tl, z). Since P(Q, z) has to be periodic: P(Q, z) = P(Q, z + d), where 
d stands for layer thickness, usually close to the full length of the constituent molecules, we 
employ the usual expansion in the Fourier series: 

p & *) = £ + £ A >< 5 « > • (4 - 7) 

l,p,q l,p,q,n 

where the summation as usual goes over < p, q < I, I > even, and n > 0, and where 

<<*£?> = j dnJ^dzP(n,z)S^(Q,z). (4.8) 
and where the new base functions are defined as 

ff£>(M = pj 1 ) ■ (4-9) 

They form an orthogonal set; 

J dhj dzS^ (fl,z)S%#(fl,z) = ^^Mm,m'**,k'*n,n' ■ (4.10) 

Similarly to the previous analysis we take (S^}} with 1 = 0,2 and n = 1 as leading order 
parameters; non-zero value of (Sqq ) and (Sq 2 ^) indicates the uniaxial smectic-A, while 
(i?2 2 ) vanishes in SmAuand becomes non-zero in biaxial smectic SmAs- In the isotropic 
phase all (S^'J) = for I, n > 0. 

In the expansion of pair direct correlation function (2.44) new terms are introduced, namely 

Otiti^z) = C 2 (f2 /-1 n) +W B Si%) +J2 W m,nS%$(ti~ 1 ti,z), (4.11) 

m,n 

where w s , w myn need to be calculated numerically, using, e.g. low-density approximation 
(3.1)-(3.2). They are integrals of c 2 (f2 &,z) with Sm,n(& &,z), as can be easily ob- 
tained with help of orthogonality relations (4.10). Since we know that the terms with angular 
momentum index I = 2 will give the leading contribution to the bifurcation, we already trun- 
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cated the above expansion accordingly; the same argument works for the Fourier series index 
n. Since the leading contribution to bifurcation comes from n — 1, we can fix it, and sim- 
plify the notation by using in further calculations Sm,n(& Q, z) = Sm}n{£l z) and 
(S$ n ) EE (S%$)- NOW we can derive the bifurcation equations corresponding to the ad- 
ditional phase transitions from isotropic, uniaxial, or biaxial nematic to uniaxial or biaxial 
smectics. However, before we present these, we show order parameters calculated using the 
self-consistent equation (2.27) with the above formulas. 



4.2.1 Order parameters in L = 2 model 



In present section we show an alternative method of obtaining the phase diagram follow- 
ing from the solutions of the self-consistent equation (2.27) in L = 2 model, i.e., we take 

(2) 

into account the expansion of C2 (4.11), calculate the coefficients Cm\ n , w s , w m ^ n by numerical 
methods (see Appendix B), and work out the non-linear integral equations for order param- 
eters. In the calculations we used low-density approximation, setting C2(f2 Cl, z — z') = 
exp[— (3 U] — 1, where U is the biaxial Gay-Berne potential (3.21). In the case of w s and w m ,n 
the calculations were significantly longer, due to the fact that the integral over z — z' has to 
be calculated in the similar way as integral over r in Cm] n , as mentioned in Appendix B. The 
equation (2.27) now becomes: 



P(n,z) = Z; 1 exp 



P 



j dti J dz'4\n l h,z-z!)P{ti,z!) 



(4.12) 



where Z s = j dfl dz exp I p J dCl dz' (f2 fl, z — z')P(Q , z' 



, and where 



cfHn 1 h,z-z') 



/ j K *m,n m,n\ J 
m,nG{0,2} 

+ w s S$(z- z')+ 



(2) -'-1" / (413) 

m,n6{0,2} 

In terms of order parameters the Eq. (4.12) can be written (similarly to Eq. (2.53)) as: 



(Sg> n ) = Z7 1 J ^ J dz *) 



exp 



/ dz'cfiti 1 (l,z-z')P(ti,z l ) 



(4.14) 
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Solving Eqs. (4.14) and Eqs. (2.52) for (A^) in the iterative manner we can calculate 
the order parameters for given density and temperature and estimate the location of phase 
transitions. 

An example of the results obtained in the way described above is presented in Fig. 4.4. It 
shows dominant order parameters for the biaxial Gay-Berne model (Sec. 3.4.2) as function of 
temperature for dimensionless density p* = 0.18 (see Eq. (3.20)). There are 9 relevant {Sm,n), 
all were calculated, but we only plot the largest ones. The range of biaxial nematic phase is 
limited by biaxial smectic, and, as can be seen by comparing Fig. 4.4(a) with Fig. 4.4(b), 
depends on the sign of energy biaxiality A e . More precisely, the region of stable N B is wider 
in case of stronger lateral interactions (where molecules are stronger attracted to their sides, 
i.e., side-to-side configuration gives the deepest minimum, see Fig. 3.5) for \ a > and A e < 
than for same signs of the biaxialities. In that respect the present approach can explain why 
Monte Carlo simulations discovered stable N B only for opposite signs of \ e and A CT [35]. 
Probably the spatially non-uniform structures gained stability earlier. As we can see, the 
transition from biaxial nematic leads directly to biaxial smectic, which is in agreement with 
[35], also L = 2 model gives ((Sqo) = (Soo)- The results obtained from the solution of the 
self-consistent equation for order parameters for I so - N u+ and N u+ -N B transitions are 
equivalent to those obtained by bifurcation analysis, the main difference being the calculation 
time, which is significantly shorter in case of Eq. (2.47) and Eq. (2.58). Thus it is useful to 
have analogical equations for the bifurcations involving SmA B and SmAjj. We present them 
in the following section. 



4.2.2 Bifurcation equations 

Given the above it is obvious that we have a larger number of possible bifurcation scenar- 
ios. Disregarding the smectic phase as a reference, we can expect five additional bifurcations; 
those are: Nb - SmA B , - SmAu, N[/ - SmA B , Iso - SmA v , and Iso - SmA B . We 
list the bifurcation equations (2.33) in the form of Eq. (2.39). They follow, as in the previous 
cases, from the expansion 




(4.15) 



which is equivalent to Eq. (2.35). The derivation is analogical to the one presented for the 
spatially uniform states. We show the equations with the help of the bifurcation matrix or 2 -*, 
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and as in case of Eq. (2.39) the bifurcation point p is determined from the roots of the char- 
acteristic polynomial, Eq. (2.41), which for the present normalizations can be written as 



det 



(2) 



2p 1 







(4.16) 




(a) Case of stronger lateral interactions, for this model Monte Carlo simulations 
showed stable N B [35] (e x ,e y , e z ) = (1.7, 1.0, 0.2) (A e = -0.06). 




r 



(b) Case of same sign of \ a and A c ; minimum in the face-to-face configuration is 
deepest, (e x , e y ,e z ) = (1.0, 1.4, 0.2) (A e = 0.04). 



Figure 4.4: Nb range for positive and negative energy biaxiality A e ; leading, dom- 
inant order parameters are plotted: uniaxial nematic (Aq 2 q), biaxial nematic (A^), 

uniaxial smectic (Sq°q), (S^q), and biaxial smectic (S^l), as calculated from L = 2 
model for biaxial Gay-Berne interaction using equation (4.14), for (a x , a y , a z ) = 
(1.4, 0.714, 3.0) (A,j = 0.58) and two values of A e . 



Preliminary results 



111 



where 1 stands for the unit matrix. From the solutions of the above equation the bifurcation 
point is chosen as the one with the lowest p , as was described previously in Sec. 2.2.2. In the 
following all the averages (Aji|„) are calculated in the reference state. 
1. Nematic uniaxial - smectic uniaxial. 

In this case the bifurcation matrix is 3x3, and the bifurcation equation reads 



/ S (0) \ 
*o,o 

b 0,0 

s {2) 

\ 0,2 / 



P0 ~ (2) 



— U 



( S (0) \ 
*0,0 

s {2) 

*0,0 

s {2) 

\ b 0,2 / 



where 



/ 



(2) 



V ^(<2> 



w ,o(A®) + wo, 2 {A$) w 2 , (A$) + w 2 , 2 {A$) 
A w 0) o + B w 0t2 A w 2 ,o + B w 2 ,2 

B W n + CqWqo #0^2,0 + C W 2 o 



and where 



35 A Q 
35 B 
35 C 



' + 10<A$) + 18<A$> 
10(Ag)+3Vl5(A ^), 

,(2)\ J_ Q/aW 



(4.17) 



(4.18) 



(4.19) 



+ 3v / 15(Ag) 
7-10(A^) + 3(A^) + 3v / 35(A^). 

2. Nematic uniaxial - smectic biaxial. 

Here the uniaxial smectic parameters Sqq, Sq 2 q, and s^l bifurcate independently, along 
Eq. (4.17), and the equation for remaining coefficients is as follows: 



( °?l > 




' 45 \ 




1=^1 














v 45 / 



where 



(4.20) 



(2) 



A)W ,o + B w , 2 A w 2fl + B w 2 , 2 

BoW n + B 2 W o B W 2 n + #2^2,2 



(4.21) 
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and where 

35j4o=7-10(A$) + 3<A$), 

70 B = 20(A<g> + vT5(A[g) , (4.22) 
70 B 2 = 14 + 20(A£>> + (A$) + V35<A$) . 

3. Nematic biaxial - smectic biaxial. 

The bifurcation eigenproblem, due to the non-trivial structure of the reference state, in- 
volves 5x5 matrix. The equation reads 



/ 


s {0) \ 

6 0,0 

s {2) 

■=■0,0 




/ 


s {0) \ 

6 0,0 

s {2) 

b 0,0 




b 0,2 


P0 r ,(2) 
— — two 

2 3 




b 0,2 




b 2,0 






b 2,0 




4,2 / 




V 


S 2,2 J 



(4.23) 



where 
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0,0 




)o,i 


(^f } )0,2 


^0 


.4 
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\ 






(-S 2) ) 


1,0 




)l,l 


(^S 2) )l,2 


-D^0,2 - Cwo,o 


Dw 2 ,2 - 






uj 3 - 






2,0 




)2,1 


(^f } )2,2 


Dw fi + Gw ,2 


Dw 2fl 4 


Gw 2i2 








Bo 




-Dw ,2 - 


C^o,o 


Dw 2 ,2 - Cw 2 fi 




iu> 2)2 - 


Hw 2 ,o 






V 


B 2 




Dwco + 


Gw 0>2 


Dw 2 fi + GU>2,2 


lu> ,0 + JU7 ,2 


Iw 2fl + 


Jw 2j2 


J 



and where 

A n = w nA^2%), B n = w s (A {2 l), 

m£{0,2} 

70C=20(Ag)-6v / 15(Ag), 
14L>=4<Ag}+3(A$>, 

70 F = -7+10(Ag)-3(A^)-3v / 35(A^), (4 . 24 ) 
70 G = 20(Ag) + Vl5<A$) + 5v / 2l(A^) , 
35 H = _ 7 + io(Ag) - 3<A$) - 3>/35<A$) , 
70 / = 20(Ag) + v / 15( A$> + 5v / 2l(Ai 4 J) , 

70 J = 14 + 20<A$) + (A$) + >/35<A$) + V35<A$) + 35(A<g) . 



Preliminary results 



113 



4. Isotropic - biaxial smectic transition. 

In this case we have the simplest bifurcation matrix, since the reference state has trivial, 
isotropic structure. The equation is as follows: 



/ s {0) \ 

b 0,0 

s {2) 

*0,0 
b 0,2 
b 2,0 

\ -S / 



6 0,0 
b 0,2 

s {2) 

b 2,0 
b 2,2 / 



(4.25) 



where 



(2) 



W„ 



V 



po 












^0,2 








W2,0 



Wo 2 















po 







w 2 ,o 



W 2 ,2 



_9_ 

Po 



7 



(4.26) 



As we can see only the order parameters (Sq 1 ^) for / = 2, 4 and n = 1 are needed to localize 
the bifurcation point. In case of - SmA B transition (2), the uniaxial and biaxial parameters 
bifurcate independently, and the equations involving uniaxial ones are the same as in case of 
N{/ - SmAu (1). It was also the case for the bifurcations where isotropic phase was taken as 
the reference for spatially uniform states (Eq. (2.46)), only this time the bifurcation equations 
involving biaxial and uniaxial smectic are different, due to the non-trivial structure of the 
reference state. This decoupling is removed for - SmAs case (3), where the bifurcation 
matrix (Eq. (4.23)) contains the elements of matrix for - SmAu from Eq. (4.17), but the 
biaxial and uniaxial bifurcations do not occur separately, the coefficients are mixed. It follows 
from the fact that the D 2 h symmetry is not broken during this transition. In case of the isotropic 
reference phase, as for spatially uniform states, the coefficients associated with pure smectic 
order parameter {Sf^), with biaxial and uniaxial smectic, bifurcate independently. 
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Preliminary results 



4.2.3 Further studies 



We have presented additional bifurcation equations for the transitions involving uniaxial 
and biaxial smectic-A phases. They extend the previously derived bifurcation equations and, 
with Eqs. (2.47), (2.58), present a way to study the stability of the biaxial nematic against 
simplest spatially non-uniform phases. The solution of those, most probably, gives similar 
diagrams, as those obtained by working out the self-consistent equation (4. 12). As an example 
we presented the temperature dependence of order parameters obtained by iteratively solving 
the self-consistent equation. The phase sequence discovered with increase of the temperature 
was SttiAb - Nb - Nu - I so. We have shown that the range of biaxial nematic phase is 
influenced by the opposite signs of shape and energy biaxiality. For stronger lateral forces 
(Fig. 4.4(a)) the temperature range 5t\ of Nb limited by SttiAb was wider than the interval 
St 2 for same signs of X e and A CT (Fig. 4.4(b)), which corresponds to the model where molecules 
are stronger attracted in the face-to-face configuration. The ratio of those ranges was found 
to be: St 2 /5ti = 0.77. In some sense it partially explains the reasons as to why the model of 
strong lateral interactions gave N B in Monte Carlo simulations [35]. This method is equivalent 
to bifurcation, but is numerically more demanding, and it is more worthwhile to consider using 
a trial one-particle distribution function analogical to (3.9), like 



trial eX P 



OfcmASkCn) + lm,nS$ n (n, 



(4.27) 



where Z tr iai ensures the normalization J dfl dz P tr i a i(£l, z) = 1, and minimize the free energy 
with respect to the parameters ot m .n and 7 m , n - This method will give more accurate phase 
diagrams. 

In the same way as we performed the analysis in search for factors that stabilize the biaxial 
nematic against I so or states, the next step would be to localize the molecular parameters 
that reduce the stability of SmArj and SmA B with respect to N B . In that way we can obtain 
a more complete picture of the region where thermotropic biaxial nematic is stable. Yet we 
already know it would not be enough. As we have mentioned, the bent-core systems exhibit a 
variety of sophisticated smectic phases, and the inclusion of those should also be considered. 



Chapter 5 
Summary and future studies 



We have presented a Density Functional Theory (DFT) study on the stability of the ther- 
motropic biaxial nematic (Nb) in comparison to isotropic liquid and uniaxial nematic phases, 
in three models of intermolecular pair potential. We started with brief description of the his- 
torical background and the research concerning biaxial nematic phase conduced so far. Then, 
we introduced the theoretical tools: the Density Functional Theory and bifurcation analysis, 
explored in the subsequent part of the thesis. We derived the bifurcation equations, in both the 
general and specific forms, for transitions involving spatially uniform states. We showed how 
the representations numbered by angular momentum index bifurcate, and that the contribution 
from pair direct correlation function to the bifurcation point comes only from the subspace of 
the order equal to the one of bifurcating representation. 

In the following part, we mainly applied DFT using the low-density approximation for pair 
direct correlation function, with the exception of the first model, where mean field theory was 
used. We focused on the determination of the phase diagram at bifurcation point for spatially 
uniform states, and commented on elastic constants and smectic phases. 

Firstly, we analysed a generic biaxial L = 2 model [10] of the pair potential. It allowed 
to acquire the exact bifurcation diagram in mean field, which emphasised how the coupling 
constants impact the transition to biaxial nematic, Landau points, and tricritical regions. Be- 
ing the simplest possible interaction giving rise to Nb, it can be used to locate the primary, 
effective factors influencing the stability of this mesophase in terms of the coupling constants. 

We continued to investigate the biaxial Gay-Berne model [36]. It was characterized by six 
constants related to shape of molecules and interaction strengths. These were combined into 
two scalar quantities describing shape and energy biaxiality, which were used to parametrize 
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the bifurcation diagrams. By varying these parameters we found the Landau(self-dual) points, 
which mark the regions where Nb is most stable in comparison to uniaxial nematic and 
isotropic liquid, and, therefore, allow to locate the microscopic parameters sets and their rel- 
ative magnitudes most important for stabilisation of biaxial nematic. We illustrated how the 
shape and energy biaxiality and the inclusion of attractive forces change Landau point posi- 
tion. The results suggest that the so-called square root rule, which relates the axes of the hard, 
biaxial ellipsoids at self-dual point, is also significant for soft biaxial Gay-Berne interaction. 
Our analysis rendered the plane parametrized by shape and energy biaxialities and defined by 
the square root rule, as a good, but not only, candidate for the approximate locations of the 
Landau points. This conclusion needs to be confirmed by more detailed studies. The biaxial 
Gay-Berne interaction was also important because Monte Carlo study predicted for it a stable 
Nb, but only for one set of potential parameters, namely for the model of strong lateral at- 
tractions, where the molecules are more attracted to their sides than faces [35]. In our studies 
N B appears stable quite easily, for many geometries, mainly because we neglect other, lower 
symmetric structures. In order to address this issue, we presented a preliminary analysis in- 
volving orthogonal spatially non-uniform phases. It indicated that the temperature range of 
biaxial nematic between biaxial smectic-A and uniaxial nematic is wider in the case identified 
in the simulations than for the model of the strongest face-to-face attractions. 

Biaxial Gay-Berne model served as the example for the way the biaxiality of molecules 
and anisotropy of intermolecular forces influence the transition to the biaxial nematic phase. 
It was also used to confirm our findings on the contribution of the L = 2 model of pair direct 
correlation function to transition points involving isotropic and nematic phases. In pursuit for 
the reasons for biaxial nematic formation in this model, we also calculated the set of bulk 
biaxial elastic constants in rod-like regime and at Landau point. Only in the vicinity of the 
self-dual point, relative values of biaxial constants changed significantly, namely the ones cor- 
responding to the deformations of one of the biaxial directors became dominant. 

A model that can be more closely related to the experimental studies where biaxial nematic 
was observed, was the system of bent-core (banana-like) molecules. We found that for banana 
with two uniaxial Gay-Berne arms the Landau point is not affected by the inclusion of attrac- 
tive forces since it was found at the same opening angle 7 as for hard bent-core molecules [38]. 
In case of three uniaxial arms, Landau point was located for 7 close to the right angle, which 
is in agreement with some of recent observations [61]. Bent-core model was also the only one 
where we introduced dipole-dipole interactions and studied their influence on bifurcation dia- 
gram. Our findings indicate that the Landau point was not only shifted towards lower opening 
angles, but also, for three parts model a self-dual line was obtained for moderate strengths of 
the dipole. In light of this, we conclude that the dipole-dipole interactions and, therefore, the 
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presence of electric dipoles in bent-core molecules are important factors in the stabilisation of 
Nb. Using the generalized Gay-Berne potential, we also took into account the model of two 
biaxial, soft ellipsoids as arms of the banana-like molecule. Bifurcation diagram showed a line 
of Landau points, indicating that the biaxiality of the arms is also important for the formation 
of biaxial nematic. 

A first step towards further studies is to investigate in detail the self-dual plane for the 
biaxial Gay-Berne model and discuss its relation to the one following from the square root 
rule. Then, one should consider the smectic structures, not only orthogonal, taken into account 
in the last chapter, but also lower symmetric, spatially non-uniform phases. In other words, 
the next(difficult) question that should be asked is: what factors enlarge the region of stability 
of biaxial nematic with respect to smectics. 

Concerning the approximations for the excess Helmholtz free energy, firstly one should 
take into consideration a more realistic reference state for the Taylor expansion, e.g., a system 
of hard molecules. Similar strategy was used in case of uniaxial molecules, and was proven 
successful in comparison with simulation results [100, 101, 124]. In those papers, a separation 
(inspired by previous studies [125] and the analogy between systems of hard spheres and hard 
ellipsoids [126]) of Gay-Berne potential into attractive and repulsive parts was performed. 
The steric part of the free energy was calculated exactly, while the contribution from the 
attractive forces was introduced in a perturbative manner. Probably, those approaches can 
be generalized to the case of biaxial ellipsoids, and some meaning of the attractive biaxial 
forces as well as a better agreement with simulations can be acquired. Also, next order terms 
in density in the Taylor expansion can be included, i.e., higher direct correlation functions 
can be taken into account. However, in this case the integrals become more challenging. 
One method of addressing this issue is a so-called weighted density approximation where the 
excess Helmholtz free energy depends on the density averaged with an appropriately chosen 
weighting function [127]. This approach was used for the case of hard biaxial molecules 
(spheroplatelets) [30, 31] and it can be extended to soft potentials. 

There are great hopes associated with the thermotropic biaxial nematic phase, concerning 
not only the promising technological applications, but also a deeper understanding of the 
mechanisms of mesophase formulation. Equally great are the experimental and theoretical 
efforts to determine the factors that increase the chances of observing this state. Present work 
is a small contribution to this ongoing search. 
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Calculation of c\ 



Here we will show that the first direct correlation function is given by 

ci Oi, [q]) = P\jjl- V ext (xi)] - In [Ag(xi)] . (A.l) 

To achieve this we will calculate the first derivative of excessive free energy functional in two 
ways, from (2.18), and then using (2.15). Differentiating (2.18) with respect to g(xi) results 
in: 

= /T 1 In [A<?(xi)] - Cl ( Xl , [g]) . (A.2) 

OQ[Xi) 

Using (2.15) we can find the following expression for the first derivative: 

6JF [g] f a _ Sn ty] Sip (x 2 ) , f j_ Sifj (x 2 ) 
Sg (xi 

where 



f , 6Q M\ Sip (x 2 ) f , 5ip (x 2 ) . . , , . 

/ dx i ttt\ x f ( + / dx * irrHefa) + ^ M , (A.3) 

J 0ip{x 2 ) dg(x 1 ) J Og(xi) 



ssi m r dx jsm_ i> M \ _ e{x2) = _ g{x2) t (A4) 



Sip(x 2 ) J \5g(x 3 ) ' J 5ip(x 2 ) 
and where we have used (2.16) to eliminate the term under the integral. We arrive at 

5T[g] 



5g(xi) 

and from (A.2) and (A. 5) we derive (A.l) 



ip (xi) = n - V ext (xi) , (A.5) 



Appendix B 
Details of the calculations 



Throughout the work we have employed certain technical methods, which were not de- 
scribed in detail, those belong in present appendix. Here we will describe the method of 
numerical integration developed to deal with the integrals that are present in Eq. (2.44) and in 
coefficients of Eq. (4.13), the methods of solving the self-consistent equations, and comment 
on the magnitude of numerical errors. 



B.l Integration procedure 

Once the bifurcation equations are formulated, we are left with only one problem: the 
estimation of the Cm]n coefficients in Eq. (2.44) (or w s , u> m ,n m Eq. (4.13)). They enter all 
of the equations for branching point (2.47), (2.58). In order to calculate them, we need a 
model for pair direct correlation function. We have chosen the low-density approximation, as 
described by Eqs. (3.1)-(3.2). Given the model two-point potential V(xi, x 2 ), we can proceed 

(2) 

to calculate the coefficients Cm,n- That task has to be performed numerically. 

We focus on the spatially uniform case; the calculations for smectic coefficients w s , w myn 
are analogical. Each Crn,n is, in fact, a six dimensional integral; 



c: 



W n = I da d(3 sm((3) / d 7 / dO I d<j> \ dr* r* 2 
o Jo Jo Jo Jo 

c(a,A 7 ,M,r*)Ag n (a,/?,7) 



(B.l) 



c(a, 13, 7 , 9, 0, r*) = exp [-V* (a, (3, 7 , 9, 0, r*)/t*} - 1 
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in the above e sets the energy unit, as in (3.20); V* = V/e and stands for the two-point 
interaction, either biaxial Gay-Berne (Sec. 3.4) or model bent-core potential (Sec. 3.5), t* is 
a dimensionless temperature defined in (3.20). In (B.l) R is a fixed maximal radius, above 
which the contribution to the integral is negligibly small. In case of Gay-Berne interaction 
the choice of R is fairly obvious, since there exists a distance at which most certainly the 
interaction is practically zero, and the above integral numerically vanishes (it follows from 
the fact that Gay-Berne potential operates on relatively short distances, and the long-range, 
attractive tail can be at some distance considered to be negligible). However in case of dipole- 
dipole interaction the choice of R is not that obvious, since the interaction is long range. In 
that case we have performed the integration (B.l) for r* < R GB , where Rgb stands for the 
radius where Gay-Berne interaction can be neglected, and for r* > R GB we have integrated 
the expanded Meyer function. So the integrand 

exp \-V GB /t* - V DD /f] - 1, for r* < R GB , 
c(a,/3, 7 ,M,0 = < 1 > (B.2) 

fV^/t* 2 , for r*>R GB . 

where V GB stands for Gay-Berne potential and Vdd f° r dipole-dipole interaction. Later the re- 
sults were verified by integrating the Meyer function symmetrized with respect to the direction 
of the dipole, where the symmetrization was imposed by the form of one-particle distribution 
function. We should keep in mind that in present approach c(a, j3, 7, 9, 0, r*) contributes via 
the Eq. (3.2): 

jj ex [P] /(N) = ~p* J dn'dnd 3 ^ Piti^n'^n, r)P(£i) . (b.3) 

This equation has physical meaning, since the equilibrium state is found as the minimum of 
T [P] = Tid [P] + Tex [P] ■ As we can see, the contribution to the above comes from sym- 
metrized c(f2 1 f2, r). The form of Eq. (B.3) is the base of all the above assumptions used in 
the calculation of c m \n- 

In (B.l) the integrand is a strongly oscillating function, and the straight forward approach 
based on the classic numerical methods is going to fail to give a precise result. All of the model 
pair-potentials V(xi,X2), studied here, as functions of r* behave as presented in Fig. 3.1(a), 
which means that c(f2 S7, r*) for fixed orientations depends on r* as in Fig. 3.1(c). We can 
see that the integral over the radius connecting centres of mass of the molecules pose the main 
problems. Therefore we have developed a "hybrid" integration method, where the 5 dimen- 
sional integral over orientations a, (3, 7, 9, and <p was calculated using Monte Carlo (MC) 
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method, and in each MC cycle, for fixed orientational variables, a precise adaptive 1 10-point 
Gauss quadrature method was employed. The number of Monte Carlo steps and the maxi- 
mum number of subdivisions of the integration interval in r* were tuned up so the relative 
error was kept to be less than 0.5%. It was estimated by calculating the integral for doubled 
number of Monte Carlo steps and interval divisions, and comparing with less precise result. 
That required a maximal number of Monte Carlo steps to be 4194304, and on average 16 

(2) 

subdivisions. In Fig. B.l an example result of the integration is presented. We show the c q 
coefficient calculated for the biaxial Gay-Berne model for biaxial, rod-like molecule as func- 
tion of the dimensionless temperature t*. Please recall that in order to solve the self-consistent 

(2) 

bifurcation equations, we need to know the values of c m \ n for given t*. 

In the case of smectic phases, the integral over z in coefficients w s , u> m ,n has to be calcu- 
lated precisely. For this reason the calculations for spatially non-uniform structures were one 
order of magnitude longer, for the integration over z was performed in the same way as over 
r* in spatially uniform case, i.e., we are left with a two dimensional integral for each Monte 
Carlo cycle. 



'This method automatically located the areas with the largest relative error, and by subsequent divisions of 
the interval in these regions and integration reduced it to the given level. 
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B.2 Self-consistent equations 

Having calculated the coefficients Cm]n for given pair potential parameters, and knowing 
their dependence on dimensionless temperature, we could proceed to solve the bifurcation 
equations. In the case of transitions from isotropic phase, the equations (2.47) are "analyti- 

(2) 

cal", in the sense that once we know Cm\ n , we can immediately find bifurcation density for 
given temperature, although not explicitly, since Cm,n do not have an analytic dependence on 
t*. In Nu-N B case the Eq. (2.58) are self-consistent equations for density depending on the 
order parameters in the reference state, which have to be calculated also in the self-consistent 
manner using (3.8). In order to perform the calculations concerning the stabilisation of the 
reference state, and the solution of bifurcation equations (2.47),(2.58) we have constructed a 
Mathematica notebook. 

The integration procedure was implemented in C language, compiled using The Portland 
Group Inc. compiler. The total CPU time consumed for the calculations, results of which are 
presented here, including the preliminary studies on smectic phases was 90000 hours. The 
numerical analysis was performed at ICM under grant G27-8, using a cluster of 1 10 nodes of 
total 292 CPUs, of which 98 machines were IBM eServer 325 (each with two AMD Opteron 
246 CPUs), and 12 were Sun v40z (each with eight AMD Opteron 875 CPUs). The work was 
financed by grant from MNiSW no N202 169 31/3455. 



Appendix C 
Symmetry adapted functions 



In current appendix we would like to present a possible derivation of orthogonal base in 
the space of real functions, which are defined on the orientation space, i.e., are mapping unit 
sphere into R. It is also a base in the space of solutions of bifurcation equations. 

In actual calculations we have used two representations of A®n functions, one with help 
of Euler angles and the other with directional cosines. Also some properties of Wigner ma- 
trices were used. We will briefly write down the formulas here. We follow the notation and 
representation of Rose [73]. 



C.l Wigner matrices 

Wigner matrices Dm,n are representations of the rotation operator in the spherical base in 
the same way the rotation matrix is in Cartesian coordinates. In present work we follow the 
convention used by M. E. Rose [73]. This section is devoted to the brief summary of formulas 
which were used to obtain the symmetry adapted version of Wigner matrices, presented in the 
next section. 

If we consider a state 4>i >m diagonalizing the square of angular momentum and one of its 
components along a given axis z, we can wonder how this state transforms under the rotation R 
about fixed direction n. From definition that transformation is described by Wigner matrices: 




(C.l) 
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where 



D (l ] ( Q)j g, 7 ) = (im' 

m m \ 117 1/ 



-iaL, Ze -if3Ly e - 



I m) 



(C2) 



and L x ,L y , L z are the angular momentum operators. 

In present study we will often encounter the relative orientation of the molecules, which 
means that we will need a formula for the Wigner matrices in that case. From definition we 
have 



R 



n 



/ m) 



(C3) 



where stands for complex conjugate of Dm,n- 

Another important property of Dm,n functions are the orthogonality relations: 

,2 



dfiD il) *(n)D (l '> ,(n 

m,n \ J m',n'\ 



(O 
m'n' 



21 + 1 



(C4) 



We will also need a rule for integrating three Wigner matrices, we can acquire it with help of 
the so-called Clebsch-Gordan series (see, e.g, [73] p. 57): 



+/x 2 ,mi+m, 2 ) 



(C5) 



where C(Zi, l 2 , I] mi, m 2 ) are Clebsch-Gordan coefficients [73], and the summation goes over 
the / satisfying the triangle rule with /i and l 2 . Now using (C.4) we can obtain the following: 



57T 



" ^v,n \ / m ,mi v / jua ;rri2 \ / 2/ + 1 



C(h, l 2 , h; Ati, H2)C{h, l 2 , W, mi, m 2 ) 



(C.6) 



C.2 Derivation of functions 

One-particle distribution function and correlation functions, in spatially uniform case, are 
defined on the orientation space. Current section is devoted to the introduction of the base in 
the space of orientation dependent, real functions, for a given model symmetry of molecules 



Symmetry adapted functions 



127 



and phase. 

Lets consider a set of irreducible tensors in spherical base. Given some symmetry 
group Q, any element of the set can be symmetrized. Namely the ^-symmetrized form of 
TJL can be written as 



(C7) 



geg 



where V (g) is the operation of symmetry element g, \Q\ is the number of elements in Q, and 
coefficients N l m ensure the normalization with respect to the following scalar product: 



T' , 

m 



1,3, 



1,3,- 



(C8) 



1,3, 



where the summation goes over Cartesian indices. T/ n can be built using the following rela- 
tion: 



m / j 
ml m2 

where Tq 1 and are as follows: 



I - I 



mi m 2 



m 



■ mi 



Ti 



m 2 ' 



(C9) 



T 1 



Ti 



± i 



±1 

7f 



(x±iy) 



(CIO) 



One of the issues that we are dealing with is establishing the relation between microscopic 
structure and macroscopic properties of the system. Usually we are given some shape of 
molecules with some model interaction, and we are left with a problem of determining the 
phase sequence. We can define the shape as invariance under some symmetry group Q B , and 
introduce a symmetry of the phase Q L . Then, we can build the sets of Q B and Q L - symmetrized 
irreducible tensors [128], and choose the base functions as following combinations: L n ■ B m 



where L, 



andB m = T' 



Gl 



Qb 



(the angular momentum index /* is fixed 1 ). If we 



assume that both Q m and Qi are equal to D 2 h , then the only non-vanishing, real, independent, 



'We refer to the "lowest, weakest or minimal coupling model" by L = 2 having in mind I* = 2 from the 
formula above. 
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non-trivial tensors remaining after symmetrization (performed along (C.7)) are: 



T 2 
x 2 



=-= (3z®z-l) , 

D 2h yfc 

1 / 

x <g> x - y ® yj 



where: 



(C.ll) 



D 2h y/2 

We arrive at the set of D 2 h -symmetry adapted functions: 
A$(A) = Lo-B =^ + ?«* (2/3) , 
A$(fi) = L B 2 = ^cos (2 7 ) sin (/3) 2 , 

A<g(fi) = L 2 B = ^cos (2 a) sin (/3) 2 , (C12) 

A^> (fi) = L 2 • B 2 = - [3 + cos (2 /?)] cos (2 a) cos (2 7) 
— cos (/?) sin (2 a) sin (2 7) , 



(C.13) 



L -B = -| + § (^3 -la) 2 , 
L -B 2 = ^|(b 1 -1 3 ) 2 - (b 2 -l 3 ) 2 

L 2 -B = ^|(b3-ii) 2 -(b 3 -i 2 ) 2 
l 2 • b 2 = (bi • ii) 2 + (b 2 ■ i 2 ) 2 - 1 (b 3 ■ i 3 ^ 2 

and where a, /3, 7 are the Euler angles [72, 73] parametrizing the rotation transforming the 
laboratory frame |i 1 ,i 2 ,i 3 | to the body frame |b 1 ,b 2 ,b 3 |. The set of symmetry adapted 
base functions (C.12) is well known [28], in present approach we can clearly see that it is valid 
for D 2 h - symmetric phase and molecules. Choosing the lowest possible I* and neglecting 
higher order terms constitutes the so-called weakest coupling model. The above allows for 
higher I* to be taken into account as well as other symmetries. We see that the first of the lower 
indices of Am in (f2) is associated with the symmetry of phase, while the other corresponds to 
the molecular symmetry group. 

Next, we will present some properties of A® n functions that are useful in calculations. 
One of them are orthogonality relations. Those are obvious when we realize that the set 
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(C.12) can be represented by Wigner matrices [73, 79] as was done in [28]: 

<T,cr'e{-l,l} 

where coefficient iV^' n ensures the normalization. In case of D 2 h symmetry /, m, and n are 
even and positive, and 



Nf=\^\ . (C.15) 



V2 



2+S a ,o+Sb,o 



Above representation is useful, since the orthogonality relations for Dm,n (C.4) are well known 



[73, 79], and we can easily see that: 



/ 



dtl Ag n (O) = N AA 5 ljV 5 m>m , 8 n>n , , N l AA = , (C.16) 



where 5ij stands for the Kronecker delta. 

We should remember that the relation (C.14), alas in present work defined only for positive 
and even I, m, n, can be easily extended to other molecular symmetry groups. Following the 
above analysis it is possible to symmetrize the spherical tensors to include symmetries lower 
than D 2 h , which for I* = 2 will mean the inclusion of other that even and/or positive indices 
m,n in Am,n- In that sense the analysis conduced using the symmetry adapted functions, 
provided they form an orthogonal base in space of real functions, can be considered general. 

The pair direct correlation function depends on relative orientation, and will be expanded 
in the above base, therefore it is of interest to take a look at Am,n fl) . From the definition 
(C.14), and using (C.3) we have 



A£„(n~n) = Jvr ,B E 

cr,cr'e{-l,l} 



(C.17) 



Using the above, and orthogonality of Wigner matrices, we can write down the following 
integration rule: 

/ dti A« n (ti'-'n) A^Un) = n aa n^' a2 iB (A) E Wn' <V • (cis) 

J a,a'e{-l,l} 



In subsequent chapters we will also need a rule for integrating three symmetry adapted base 
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functions which follows from relation (C.6): 



J ^A» jni (fi)Ag jn2 (fi)AS in3 (fi) = Al_ N ^n 1NT ,n 2NT ,n 3 

y^i,m 3 ,ms(M), I G {2,4} 



(C19) 



where we have denoted 



iVmi,m2,m 3 ({crj}) = (— l) 
11,112,13 



^(T3 7772+05 7713,-0"! 7771 ^ (TA 772+<X6 773,— CT 2 771 

C (2, 2, Z; o- 3 m 2 , <7 5 m 3 ) C (2, 2, Z; a 4 n 2 , <7 6 n 3 ) 
ze {1,2,3,4,5,6}, 



(C.20) 



and where C (li, l 2 , 1; m 1 ,m 2 ) are Clebsch-Gordan coefficients [73] in decomposition of state 
I, m on to li, mi an d l 2 , ^2, an d the summation goes over all the as, each taking values of 1 
and —1. 

There are many ways of introducing the base functions. Most of them take into account 
the analysis of the symmetries involved in the problem. We have presented one possible 
derivation which stresses out the way in which phase and molecular symmetry come into the 
picture. The other way would be to postulate (C.14) by requiring the proper behaviour under 
the D 2 h symmetry and realness, as was done for instance in [28]. 

Main benefits from using the irreducible representations in present work come from clear 
transformation laws under rotations, built-in orthogonality, and completeness. The Cartesian 
tensors usually come in reducible form, and therefore present a more difficult objects when it 
comes to actual calcualtions. On the other hand, since we are more used to the Cartesian form, 
some interpretation and meaning can be retrieved more easily in that base. 



Appendix D 



One-particle distribution function and 

order parameters 

During phase transition a symmetry change takes place. In isotropic - uniaxial nematic 
case the rotational invariance of the disordered liquid is broken, and replaced by the symmetry 
under the action of elements of group containing the rotation around the director and the re- 
flection in the perpendicular plane. One possible method to model that type of behaviour is to 
define some quantities that in higher symmetric phase vanish, while in more ordered structure 
obtain a non-zero value. Those are usually defined as averages over an ensemble of the base 
functions (C.12), and are called order parameters. 

Of course the symmetry change happening over a phase transition has to influence the one- 
particle distribution function. If we employ the order parameter approach it seems natural that 
the distribution should be parametrized by those objects. We will begin with an example of 
such approach. 

Lets take into account two orthogonal tripods, one associated with laboratory, usually 
chosen to coincide with the director frame: jii,i 2 ,i 3 j, and the other with the molecule: 

|bi, b 2 , b 3 |. For the moment we take into account the uniaxial nematic phase, and molecules 
having the symmetry of ellipsoids of revolution, and write the one-particle distribution func- 
tion as: 

P(\l u h l })=6 (b 1 -i 3 ) 2 + r 7 (b 2 -l 3 ) 2 + tf (b 3 -i 3 ) 2 . (D.l) 
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One-particle distribution function and order parameters 



The symmetry requires that in the laboratory frame there is a rotational invariance around the 
director, and that the distribution is not changed by the reflection in the perpendicular plane. 
That is the reason why only 1 3 is present in the above formula, and that the expression is 
quadratic in products of the base vectors. We can find a transformation to the orthogonal base, 
and receive 



PiUub* 



1 3 f x 2 



+ b 



bo • I 



(D.2) 



Now the parameters a and b can be calculated as ensemble averages of the expressions in 
square brackets, using the orthogonality of the new base functions, and identified with the 
order parameters. The expansion of the form (D.l) is inspired by the fact that one-particle 
distribution function undergoes a qualitative transformation on line leading from one state to 
another of lower symmetry. It is of course a very simple approximation. There are several 
others available. One of them, clearly inspired by the form of the self-consistent equation 
(2.27) is to assume 



P(n)~exp (X^n^nCA) 



(D.3) 



\l,m,n 



Another is a Gaussian model, where one-particle distribution function is taken to be of the 
following form: 



P(a, (3, 7) ~ exp [a (a 2 + (3 2 + 7 2 )] . 



(D.4) 



Expression in Eq. (D.l) is the simplest possible model of one-particle distribution function, 
and should be treated with care. One obvious flow is the fact that probability of this form is 
not always positively defined, and may become negative. Despite that problem (and others) 
the approximation works surprisingly well. The Gaussian (D.4) and exponential (D.3) models 
provide a better behaviour of P(fl) but the relation of coefficients a and cn m n to the order 
parameters is not as simple as in case of Eq. (D.2). 

Having performed the orthogonalization that led to the expansion (D.2), we have changed 
the base to a certain class of functions which by construction are invariant under the operation 
of Dooh symmetry. Using some general assumptions concerning the phase and molecules, we 
have received the uniaxial version of the expansion in the symmetry adapted base introduced in 
(C.12). In case of the molecules orientations of which, due to symmetry, need to be described 
by three vectors, the same method can be used to obtain the whole base set in D 2 h - symmetric 
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case. It should be noted that in the initial expansion (D. 1) the higher order terms in directional 
cosines, which in the language of Am,n would mean the higher I, are not a priori forbidden. 
The choice of the lowest order terms is a part of the approximation. 

It is natural to choose the order parameters as averages of Am,n functions, we can write 

P W = ^2 + E ^ <A8n>A£„(n) • (D.5) 

l,m,n 

Using the relations in Eq. (C.16) we can immediately see that 

(A«„} = J dflP(n)A^ n (Q), (D.6) 

where as usual the average of quantity A is defined as (A) = J dfl A(Cl)P(Cl). From the 
above set (Aq 2 q) is a well known uniaxial order parameter present in Maier-Saupe theory [6]. 
Second one, (Aq 2 ^), was introduced by Freiser [5], (A^q) was seen in the paper by Alben, 
McColl, and Shih [81], and the last but not least, (A^) was for the first time defined by 
Straley [10]. 

The equation (D.6) can be treated as the definition of order parameters, once a model of 
one-particle distribution function has been chosen and the equilibrium state found. 
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Mean field approximation 



There exists a class of models that treat systems of interacting molecules by certain ap- 
proximations concerning the potential energy and distribution functions. Common to all of 
those is the presence of a pair interaction, one-particle distribution function, and an effective, 
average potential that models all of the forces that a given molecule is subjected to. They are 
called mean field theories. Presently, we show a way in which they are introduced. 

As was mentioned in Chapter 2, the mean field approximation can be introduced as a 
further expansion of two-point correlation function modelled by c(x±, £2) (see Eq. (3.2) and 
(3.1)). In current appendix, we will present a different method of derivation. 

Let's consider a system of N molecules interacting by N -particle potential V (x±, . . . , x^) 
and an arbitrary iV-particle distribution function P^ = P (xi, . . . , xn)- We can write down 
the free energy as a functional of Pjv, namely 




(E.l) 



If we assume that the iV-particle distribution function can be taken as a product of one-particle 
distribution functions (in other words we assume the independence of distributions, which is 
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Mean field approximation 



not necessarily true), the above changes into: 



F[P]/N=- I dx 1 



dx 2 P(xi)V (xi, x 2 ) P(x 2 ) 



f k B T I dx 1 P{x 1 )hiP{x 1 ) . 

(E.2) 



Now we can perform the minimization of the above, with normalization condition for P(xi), 
to acquire the following equation for equilibrium distribution 



Peg Ol) = Z ^Xp 



-(3 J dx 2 V(x 2 ,x 1 )P(x 2 ) 



(E.3) 



where V(x 2 ,Xi) is the pair potential, and where Z is a normalization constant so 
J dx\ P eq (xi) = 1. As can be seen from the above equation, given molecule is subjected 
to an average effective potential V e jf (xi) = J dx 2 V (x 2 , x x ) P(x 2 ). 

Eq. (E.3) is well known, and widely used. Its success lies mainly in the relative simplic- 
ity, compared to the variety of systems and behaviour that can be modelled by it. However, 
clearly, the competition of the entropic and potential energy terms in Eq. (E.l) is disturbed in 
Eq. (E.2), which results in overestimated transition temperatures obtainable from Eq. (E.3). 
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